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INTRODUCTION 

The irreversible annihilation reaction A+A — > 0, also known as the mutually annihilating 
random walk, is a fundamental model of non-equilibrium physics. The particles A perform 
chaotic motion due to diffusion and after the mutual collision they may react with constant 
microscopic probability K per unit time. Usually it is assumed that resulting molecule 
is inert, i.e. chemically inactive and without any backward influence on the motion of 
reacting A particles. Many reactions of this type are observed in diverse chemical, biological 
or physical systems. For instance various models such as formation of domain in some 
magnetic materials [1], annihilation of excitons in crystals [2] or model for spreading of 
opinion of voters in one dimension [3] can be described in terms of annihilation process of 
this type. 

The usual approach to the problems dealing with chemical reactions is based on the use 
of the kinetic rate equation [4] for concentration field n(t,x). It leads to a self-consistent 
description analogous to the mean-field approximation in the theory of critical phenomena in 
the sense that fluctuations in the concentration are neglected. Equivalently one can assume 
that the particle concentration is spatially homogeneous n = n(t). This homogeneity can 
be thought as a consequence of either very high mobility of the reactants or of a very small 
probability that a reaction actually occurs after mutual collision of reacting particles. For 
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the annihilation process A + A — >■ kinetic rate equation can be formulated as 



d t n(t) = -K n 2 (t), 



(1) 



which can be easily integrated and the obtained solution is 



n(t) 



1 + n K t ' 



(2) 



where no = n(0) is initial number of particles. For long-time (t — > oo) asymptotic decay 
equation (1) predicts power law behaviour for concentration n(t) ~ t^ 1 without any depen- 
dence on the value of space dimension. This is a common situation observed in the mean 
field-like theories. In what follows we will refer to the value of time exponent a in power 
law dependence for concentration n(t) oc t~ a as the decay exponent. Note also that the 
long time behaviour does not depend on the initial number no of reacting particles. In the 
other case, when the particle mobility is sufficiently small, or equivalently, if the microscopic 
reaction probability Kq becomes large enough (so particles react immediately after mutual 
collision) there is a possible transition to a new regime. In it is more probable that the given 
particle reacts with particles in its neighbourhood than with distant ones. This behaviour 
is known as the diffusion-controlled regime [4, 5]. To gain physical insight let us consider 
diffusion process (also known as continuous random walk) to be responsible for the motion 
of particles. A well-known property of diffusion [6] is the re-entrancy of the visited sites in 
low space dimensions. In particular, for d = 1 and d = 2 the probability that the diffusing 
particle will ever return (t — > oo) to the starting point is equal to 1. Physically it means 
that the diffusing particle sweeps thoroughly its local neighbourhood and thus it is highly 
probable that it will react with another particle in its vicinity. Hence, it is reasonable to 
expect that after short period of time the system will be in a state, where there is a lot of 
isolated particles, that need effectively longer time to traverse to each other and hence to 
annihilate. This mechanism can effectively slow down the time evolution of the process and 
thus lower the decay exponent to other value than 1 predicted by equation (2). The approx- 
imate value of the decay exponent can be guessed according to following scaling argument. 
The re-entrancy property leads to the scaling relation 



where ~ denotes corresponding scaling relation between physical quantities. The root mean 
square distance for the diffusing particle scales as r(t) ~ (Dt) 1 ^ 2 and therefore the mean 



V(t) ~ r d (t) , 



(3) 
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particle number should behave as 



n(t) oc 



1 1 1 



(4) 



V(t) t d / 2 t!+ A ' 



where the exponent A denotes the deviation from the space dimension 2 via relation 



For the space dimension d = 3 we have V(t) ~ t, because now the diffusing particle effectively 
explore always new volume and the re-entrancy property can be neglected. Therefore the 
same behaviour as the one described by (f ) would be observed. >From this simple analysis 
it could be estimated that space dimension d c = 2 is the upper critical dimension for the 
annihilation process, above which the mean field approximation is valid. A more rigorous 
treatment [7] based on renormalization group proves this conclusion and also produces a 
logarithmic correction for n(t) at the critical dimension, which could not be determined by 
the simple scaling analysis. In the preceding discussion we have considered only the diffusive 
motion of reacting particles. However, a typical reaction usually occurs in liquid or gaseous 
environment. Thermal fluctuations of this environment or some external advection field such 
as atmospheric eddies could have additional influence on motion of the reacting particles. 
Therefore, it is interesting to study what effect the external velocity field can have on the 
annihilation process. 

The most flexible approach to the theoretical analysis of the effects of fluctuations in 
reaction kinetics seems to be the second-quantization method due to Doi [8]. Most of the 
renormalization- group studies of the effect of random drift on the annihilation reaction A + 
A — > in the framework of the Doi approach have been carried out for the case of a quenched 
random drift field. Potential random drift with long-range [9, 10] and short-range correlations 
[11] have been studied as well as "turbulent" flow (i.e. quenched solenoidal random field) 
with potential disorder [12, 13]. For a more realistic description of a turbulent flow time- 
dependent velocity field would be more appropriate. In Ref. [13] dynamic disorder with 
a given Gaussian distribution has been considered, whereas the most ambitious approach 
on the basis of a velocity field generated by the stochastic Navier-Stokes equation has been 
introduced here by two of the present authors [15]. >From the point of the Navier-Stokes 
equation the situation near the critical dimension d c = 2 of the pure reaction model is even 
more intriguing due to the properties of the Navier-Stokes equation. It is well-known fact 



d = 2 + 2A . 



(5) 
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[16] that in the case of space dimension d = 2, there is inviscid conservation law of enstrophy 
absent in the three dimensional case. Calculations in Ref. [15] were performed in the one-loop 
approximation. As may readily seen from examination of the Feynman graphs, in the one- 
loop approximation there is no influence of the velocity fluctuations on the renormalization 
of the interaction vertices. However, the influence of higher order terms of the perturbation 
series can have significant effect on the critical properties. 

To this end, the most suitable choice for generation of random velocity field is the stochas- 
tic Navier-Stokes equation, which can be used to produce a velocity field corresponding to 
thermal fluctuations [17] and a turbulent velocity field with the Kolmogorov scaling behaviour 
[18]. 

A powerful tool for analyzing asymptotic behaviour of stochastic systems is provided by 
the renormalization-group (RG) method. It allows to determine long-time - or infra-red 
(IR) - asymptotic regimes of the system and also it is very a efficient tool for calculation 
of various universal physical quantities, e.g. critical exponents. The aim of this study is 
to examine the IR behaviour of the annihilation process under the influence of advecting 
velocity fluctuations and to determine its stability. Using the mapping procedure based on 
the Doi formalism [8] an effective field-theoretic model for the annihilation process will be 
described in detail in the Sec. 1. Consequently, the RG method is applied to this model 
and within the two-parameter expansion the renormalization constants and fixed points of 
the renormalization group are determined in the two-loop approximation. The non-linear 
integro-differential equation, which includes first non-trivial corrections to the (1), is obtained 
for the mean particle number and it is shown how the information about IR asymptotics can 
be extracted from it. 

Another aspect of the annihilation problem is connected with its theoretical description 
in the form of a functional integral with a given action. This action resembles actions for 
field-theoretic models of critical dynamics obtained from the Langevin equation [19] within 
the Martin-Siggia-Rose approach [20]. 

In the Langevin equation the random field compensates for dissipative and reactive losses, 
hence bringing about a steady state of dynamics of the system. In reaction kinetics the ran- 
dom sources and sinks, in fact, reflect the real physical situation, in which during the chemical 
reaction of follow-up species, particles can appear or disappear due to uncontrolled random 
interaction with a particle bath, e.g. due to active chemical radicals. The interpretation of 
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random fields in the Langevin equation as physical sources and sinks is rather problematic. 
Therefore, we propose to analyze the alternative approach provided by the master equation 
with terms corresponding to interactions with the bath. 

In the present paper we follow the ideas of [4] and describe the random sources and 
sinks in terms of new birth and death reactions in the master equation for the single-species 
annihilation reaction. The simplest choice does not conserve the particle number and we 
have no possibility to compare it with the standard Langevin approach. A slightly more 
involved set of birth-death reactions allows to conserve the particle number and the result 
may be compared with that of the standard multiplicative noise in the Langevin equation. 

The aim of this paper is to give a detailed account of how the field-theoretic approach 
can be applied to the analysis of the large-scale asymptotic behaviour of annihilation process 
A + A — )> 0. We give an elaborate derivation of the action functional for this reaction and 
show how the velocity fluctuations and effects of sinks and sources can be described. For the 
analysis of the large scale behaviour methods of the renormalization group are employed. 

This paper is organized as follows. In Sec. 1 a detailed description of the Doi approach 
for construction of the field-theoretic model for the annihilation process A + A — > is 
summarized. In Sec. 2 it is shown how both the Kolmogorov scaling and thermal fluctuations 
can be included into the model and the ultraviolet (UV) renormalization of the model and 
elaborated algorithm for the calculation of the renormalization constants is described. Fixed 
points of the RG are classified together with their stability regions and possible scaling 
regimes are presented and the integro-differential equation for the mean particle number is 
derived and analysis of its solution is given. Sec. 3 is devoted to the introduction of random 
sources and sinks into the field-theoretic model of reaction-diffusion processes. We recall the 
basic features of the Doi formalism and construct the basic field-theoretic dynamic action 
functional together with scaling analysis of the dynamic actions is performed. 

1. FIELD-THEORETIC REPRESENTATION OF THE MASTER EQUATION 

1.1. INTRODUCTION 

In our work we are mainly interested in the annihilation process A + A 0, therefore 
in this section we describe a derivation of the field-theoretic model for such process, which 
allows to take into account spatial inhomogenities and randomness in individual reaction 
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events. Probably the most fundamental description of reaction processes is based on the use 
of master equation [4]. We summarize general principles of how such master equation can be 
mapped onto a suitable Fock space, which makes it possible to use very powerful methods 
of quantum field theory. The resulting action functional can be treated systematically by 
such methods as Feynman diagrammatic technique, renormalization group and operator 
product expansion [21]. Here we focus mainly on the field-theoretic description of diffusion 
and reaction process, but it is possible [22] to generalize this approach to include effects as 
multiple reaction schemes, disorder effects, influence of spatial boundaries etc. 

Let us start with the particles on a regular, infinite, hypercubic lattice with the lattice 
spacing a in d-dimensional space. The sites of lattice can be labeled by natural numbers 
(i = 1,2,...). The A particles are performing continuous random walk on this lattice 
(random hopping betweeen adjoint sites) with the diffusion constant D/a 2 (the factor a 2 
will be eliminated in the continuum limit). It is assumed that reaction process can happen 
only for particles that occupy the same lattice site with the probability rate K . A complete 
(microscopic) description of such a stochastic problem can be given in terms of evolution 
equations for probabilities P(t; {n}), where {n} is given microstate {ni, ri2, ■ ■ ■} characterized 
by rii particles at site 1 , ri2 at site 2 and so on. These equations (known as master equations) 
express the balance between incoming and outcoming probabilities [4] for the state {n} and 
in a compact notation they can be written as 

dP{t ' d l n}) = E R m^nP(t; {n}) - £ R n ^ m P(t; {n}), (6) 

{m} {m} 

where R m -+ n is the transition probability rate from the state m to the state n. According 
to the work of Doi [8] (see also [23]) such a system of coupled differential equations (6) 
can be rewritten in terms of creation and annihilation operators well known from quantum 
mechanics. If there is no site occupation restriction, bosonic operators for each lattice site i 
can be introduced with the following commutation relations 

[Si, at] = 5ij, [di, dj] = [at, at] = 0. (7) 

The ground state |0) is defined as 

Oj|0) =0 for all sites i, (8) 

which corresponds to the empty lattice (without any A particle). >From the bosonic com- 
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mutation relations (7) important relations follow: 

a?aj = nar 1 + a\a n , a^f = naf _1)t + a nt a, . (9) 

The state |{n}) with the given lattice configuration {n} = {ni,ri2, . . .} is introduced with a 
normalization different from that used in the second quantization method in quantum field 
theory 

|{n»=a^...|0>. (10) 
Using relations (9) it can be directly shown that 

CLi\{n}) = n i \{n 1 ,n 2 , . . . - 1, . . .}), 
«!lW) = \{n 1 ,n 2 ,...,n i + l,...}), 

a\ai\{n}) = rii\{n}) , (11) 



where the last relation legitimizes the identification 

At A 



(12) 



for the number operator at site i. The scalar product between two states \{n}) and \{m}) 
can be obtained 

i=\ 

where 5{j stands for Kronecker symbol and n\ = 1 x 2 x . . . n is the factorial function. The 
complete information about the stochastic system is embodied in the probabilities P(t; {n}) 
and in the Doi formalism it is incorporated into the state vector |$), which is defined as 
follows 

|<&(f)> = E P ^' W)l W> = E P ^' M) a\ ni a ] 2 2 ■ ■ ■ |0>, (14) 
{»} {«} 

where the sum runs over all possible lattice occupations. Now the task is to rewrite master 
equation (6) into the Schrodinger-like form for the state vector |<E>) 

!l$(f)> = -H\*{t)h (15) 

with some "Hamiltonian" H, whose exact form depends on the system under consideration. 
Then the equation (15) can be formally integrated to obtain = e~ H '|$(0)). The initial 

state |$(0)) has to be specified for the full description. In the case of chemical reactions 
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initial distribution of particles P(0; {n}) is usually prescribed and initial state follows from 
the definition (14). From the technical point of view the most convenient choice for the singe- 
species annihilation reaction is the Poisson distribution. Unlike the bimolecular reaction [24], 
the long-time behaviour is independent of the concrete form of initial conditions. 

Let us derive the diffusion part Ho of the Hamiltonian H, which corresponds to the 
diffusive (random walk) movement of A particles. First consider a two-site system with n\ 
particles at site 1 and n 2 particles at site 2 with the one-directional hopping process 1 — > 2 
at the rate D /a 2 . For such a process the master equation (6) is 

dP{n 1 ,n 2 ) D . Dq . 

= — (m + l)P(ni + 1, n 2 - 1) -7iiP(ni, n 2 ) , (16) 

at a z cr 

where ni + 1 and n\ are combinatorial factors resulting from the fact that A particles jump 

independently of each other. Multiplying both sides of the equation (16) by the term a^a^ 1 2 

and performing sum Yl{m n 2 } over an P oss ible occupations of sites we arrive at 

^ = ^(ajai-ata 1 )|$>. (it) 

at a z 

This result is easily generalized to the two-directional case 1 -H- 2 described by the following 
master equation 

^ = -^(4-at)(a 2 -a 1 )|$) (18) 
at a z 

and hence the diffusion part H D between two given sites can be written as 

H D = ^(di-d\)(d 2 -d 1 ). (19) 

(X 

Now consider the annihilation process at given site 1 and derive the corresponding part Hr 
of the Hamiltonian. Because any two particles can react together, we can write 

= K {n + 2)(n + l)P(n + 2) - K n{n - l)P(n), (20) 

where again the combinatorial factors are taken into account. After short algebraic manip- 
ulations and the use of the relations (7) and (11), relations (20) can be rewritten in the Doi 
formalism as 

^ = K o (a 2 -a t2 a 2 )|0), (21) 
from which we deduce the reaction part of the Hamiltonian 

H R = -K (a 2 - a t2 a 2 ) . (22) 
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Results (19) and (22) are easily generalized to include all the lattice sites. Consequently, 
the total Hamiltonian that accounts for diffusion and annihilation process on the hypercubic 
lattice has the following form 

Hd + Hr = ^J2 $ - a 5) & - ">) - K ° - *?® » ( 23 ) 

<ij> i 

where the first sum runs over the neighbouring sites % and j. Non-Hermitian Hamiltonians 
such as one given in (23) are often observed in the case of systems out of equilibrium [3], which 
cannot be obtained as dynamical counterparts of some static models. Non-hermiticity also 
means that the reaction rates in (6) do not satisfy the detailed balance condition and thus 
the equilibrium state cannot be characterized by the Gibbs distribution. The Doi formalism 
also exhibits other differences from the usual quantum mechanics. They are caused by the 
fact that physical observables cannot be given as bilinear products since according 

to (14) this would imply expressions bilinear in probability P(t; {n}). 

Let us now take a closer look at derivation of the ensemble average value for an observable 
quantity A within the Doi approach. It is physically reasonable to assume that A can be 
expressed as a function of the occupation numbers A = A({n}). Examples of such quantities 
interesting for the case of chemical reactions are 

(a) mean particle number (concentration) 

n ^Y,a\ai (24) 

i 

(b) two-point correlation function (between sites i and j) 

C(i,j)< — > alaiOjCLj (25) 

The ensemble average of A is then clearly given by the expression 

(A(t)) = Y,P(t;{n})A({n}) (26) 
{«} 

and from the technical point of view it would be very convenient to have a projection state 
{V\ such that following identity is valid 

(A(t)) = P (t; {n})(V\A({^a})a^ap . . . |0> = <P|i|<l>(f)>. (27) 
M 
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Substitution of the formal solution of the "Schrodinger" equation (15) leads to yet another 
form of the last expression on the right side of (27) 



{A(t)) = {P\Aexp(-Ht)\*(0)). 



(28) 



Here the operator A({a)a}) is obtained from the classical function A({n}) with the substi- 
tution rii — > a\ai at every site i, what is justified by (12). Hence by comparing (26) and (27) 
it is easy to guess that the projection state (V\ should satisfy relation 



(V\A({^a})a\^a^ . . . |0> = A({n})(V\^S^ ■ ■ ■ |0> = A({n}) 



It implies that the following two conditions 



(29) 



(p\al = (P\ for every site i, (V\0) = 1 



(30) 



have to be valid for the operator V. By direct use of relations (7) we can conclude that the 
following choice of the projection state 



{V\ 



exp 



E' 



(31) 



serves our purpose. Note that if the operator A is written in the normal-ordered form 
(all creation operators are commuted to the left), it can then be written in terms of the 
annihilation operators only using properties (30) of the projection state (V\, i.e. 



{V\A{{a)a}) = (V\N \A N ({a\a})} = (V\A N ({l,a}) 



An^cl^cl}) 



(32) 



where iV 



The normal form An of the operator A is defined by relation A = N 
denotes the normal product (creation operators are put to the left of annihilation operators). 
The a-dependent operator A N ({l,a}) corresponding to (24) is then simply hi, while the 
correlation function (25) corresponds to operator ai&j + SijCij. For technical reasons it is 
reasonable to commute the factor e^ ai to the right in equation (27). In order to do this we 
employ the following formula 

e^a\ = {a\ + l)e a % (33) 
which is derived from relations (9) and thus equation (28) can be cast into the form 



(A(t)) 



A N ({l,a})exp(-H({tf + l,a})t) 



e E ^ l $(0)), 



(34) 
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where the substitution a\ — > a\ + 1 was performed both in the expression for the operator A 
(with the subsequent substitution a\ — > 0) and also in the original Hamiltonian H{{a\a}). 
For further use let us write expression for the mean particle number (24) 



n(t) = ( 



^2 o-i exp(-F({a f + 1, a})t) 



(35) 



In the case of annihilation process described by the equation (23) thus we obtain 

<ij> i 



(36) 



As was mentioned above, the convenient choice for the initial condition is the Poisson dis- 
tribution, which for a given site i corresponds to 



p(m) = e 



rii 

-n n 

nA 



(37) 



where n stands for the mean particle number. Using definition (14) we arrive at the initial 
state vector in the form 



|$(0)> = £V no ^ T a t 1 ni Y,*~ n —^\ n2 ■ ■ ■ 1°) = n e " n ° en0EiaI l°)- ( 38 ) 

{ni} {n 2 } i 

We see that after the substitution a) — > a) + 1 the term e~ n ° drops out. 



1.2. CONTINUUM LIMIT 

In the field of critical phenomena emphasis often lies in the analysis and determination 
of possible behaviour of the studied system. It turns out that near a second-order phase 
transition [25] a whole set of different models behave in the same way. Despite the fact that 
they describe different physical systems, they can exhibit the same behaviour of so called 
universal quantities. They are model-independent, but can depend on universal parameters 
such as space dimension, number of components of the order parameter or symmetries of the 
system. Examples of such universal quantities are critical exponents [26], describing singular 
behaviour of various functions. 

Now we summarize the main points of the derivation of the continuum limit for the 
Hamiltonian (23), which allows us to study universal properties of the annihilation process 
A + A^-0 around its critical dimension. In fact, here we just briefly outline a few important 
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steps following [22] of the introduction of the functional integral of the continuum limit by 
the popular interpolation procedure. The alternative operator approach will be sketched 
in Sec. 2.2.1. The main task consists of evaluation of matrix elements for the evolution 
operator exp(— H({a\ a])t). In order to do this we apply the Trotter formula [27], according 
to which the exponential e~ Ht can be written as the infinite product 

exp(-Ht) = km (1 - HAt) t/At = (1 - HAt)(l - HAt) . . . (39) 

Here we assume that NAt = t and at the end of our derivation we let the number of time 
slices N — > oo (or equivalently At — > 0). Now into each time slice complete set of coherent 
states is inserted, which results into mapping of operators a], di onto complex numbers. 
Coherent states are explicitly defined as [27] 

|^>=exp(-i|vf + ^|0> (40) 

and they form the eigenstate basis of the annihilation operator 

a|V> =#/>>, (^=^1, (41) 

where the star stands for complex conjugation. The important property of coherent states 
is their overlap function between different eigenstates 

<^ 2 > = expf-i|Vi| 2 - \m 2 + r^ 2 ) ■ (42) 



T T 1 2 

For a single site we can write an identity resolution in the form 

l = E i [^H = E i [l n )H^= / —1^1, (43) 

n m,n 

where the orthogonality relation 

J dip* dip exp(-\ip\ 2 )ip* m ip n (44) 



ixm\ 



has been used (note the apperance of the weight function e '^' 2 , the measure adopted here 
is dip* dip = Re^Im^). Generalization of (43) to the whole lattice is straightforward 



n 



diP* {ij) diP {id) 



TT 



(45) 
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where now {ip}j = (i>(i,j),ip(2,j), ■ ■ •) denotes the set of all eigenvalues corresponding to the 
annihilation operators a, at each lattice site at the time instant jAt (j — 0, . . . , N). Inserting 
(45) into each time slice in (39) we get 



exp(— Ht) = lim 



N 



Wiv)(n(w. 

•3 = 1 



exp(-H({a\a})At) 



where Af is the normalization constant and we have introduced the notation 



Mo 
(46) 

(47) 



for the functional measure. Note that if we deal with normal ordered Hamiltonian (which 
we shall assume and which explicitly is the case for (36)), then using relations (41) we can 
immediately write 



exp(-H({a\a})At) 



j'-i 



)eM-H({rh: Wi-i)At), (48) 



where is obtained by the replacement of operators by their eigenvalues 

a-i — > oj — > i>i- The remaining term in (48) more precisely stands for the expression 



3-1 



n<*< 



(ij) 



^(ij-l) 



(49) 



which, with the of the overlap relation (42), can be rewritten as 



VW-i)^= ex P^-^(* j)tV'(ij) - ^j-i)]^expQ| 



(ij-l)l 



(50) 



The whole scalar product from (46) can be expressed to the first order in the time increment 

as 

n(^<«>^W-i))=exp(^ (51) 

so that in the continuum time limit At — > we obtain the functional integral representation 
for the evolution operator in the form 



exp( 



-m = / 



Af 



Wn)x 



expQ|^(0| 2 - ^l^(0)| 2 - f dtWMi + H({^})]j (Mo 



(52) 
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In calculation of expectation values such as (24) and (25) we have to act on the expression 
(52) from the left by the projection state (V\ and from the right by the initial state |$(0)). 
It can be done in the following manner: first, we note that the following relation holds 



$ (0)\nexpf-i|^(0)| 2 )= exp(5> ^(0) - n - |^(0) 



(53) 



where we have used the initial state in the form (38) and relations (40,41). Using equation 
(33) we proceed as follows 



V 



e£iaie£i(-*lifc(t)l a +ik(t)aj) 



1 



JJexp(oi)exp(V>i(*)a£ 



i k=0 



exp(aj)a. 



tfc 



Y[ exp I ipi(t)(al + 1) j exp(a. 



(54) 



Putting together terms from (52-54) and inserting them into expression (28) for the expec- 
tation value of the quantity A we arrive at the important expression 



(A(t))=M- 1 



A N ({l^})exp[S({r,m, 



(55) 



where the functional A N ({l,tp}) is obtained from the normal form A N ({1, a}) of the operator 
A by replacing the operators a; by their eigenvalues ipi and the action functional S is given 



as 



S({r,ip}) = +"o#(°) - no - |Vi(0) I 2 - jf 



+ ■ (56) 



The normalization constant A/" is now fixed by the condition 



(57) 



Before taking the continuum limit in space let us note that the initial term — |-^j(0)| 2 in 
(56) can actually be dropped in calculations within perturbation theory. In fact, we will 
define the functional integral in (55) in terms of perturbation expansion and use the bilinear 
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part of the dynamic action to generate propagators, since otherwise the convergence of the 
functional integral is somewhat problematic [22]. Thus, we arrive at perturbation theory with 
a retarded propagator, which we choose such that its value at coinciding time arguments is 
zero by definition. In the traditional field-theoretic parlance this is tantamount to defining 
the time-ordered product at coinciding time arguments as the normal-ordered product. Since 
we are dealing with expectation values of functions of the annihilation operators only, a short 
reflection of the perturbation expansion reveals that all graphs containing vertices brought 
about by the initial term — 1^(0) | 2 either are proportional to the equal-time value of the 
propagator, which we have chosen to vanish, or contain closed loops of propagators and 
therefore vanish as well. 

Continuum limit then can be performed in traditional manner according to the substitu- 
tion 

-> ^ ^(*> x ) fld > W (*) ^ f (*> x )' n o ->■ n a d , (58) 

(X 

i 

that leads to the field-theoretic action for the scalar fields ^(x, £) and t/;(x.,t) 



S = - 



J^dt J rfx|^ t ^- J D ^ t V 2 ^-A [l-^ t2 ]^ 2 }+ J dx 



ip(t, x) +n V ;t (0,x) -n 

(59) 

It corresponds to the lattice Hamiltonian (23). After the shift ^ -)• ^ + 1 we arrive at 
the desired result, which is the the field theoretic action S for the annihilation reaction 
A + A ->■ 



S = 



(60) 



which corresponds to the continuum limit of the action used in (36). It should be noted 
that the shift ^ — > ^ + 1 allows to replace the final term J dxi/)(t,x) in (59) by the initial 
term J cfoc0(O,x) and the contribution of the latter to the perturbation expansion vanishes 
by the same token as the contribution of the quadratic initial term. Therefore, we have not 
included this linear initial term in the dynamic action either. 

The action functional (60) will be used in calculation of physical quantities such as the 
mean value (34) that can be expressed as the functional integral 

(A(t)) = J V^V^A N {1, ij(t)} e s . (61) 
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Here Vip^Vip is the continuum functional measure. The continuum limit of the expression 
(35) for mean particle number becomes 



ra(t,x) = (0 



V>(x) exp(^-H{^(t) + l,ip(t)}?jexp(no J dx ip ] 



0). (62) 



J. 3. L ANGEVIN EQ UATION AND FOKKER-PLANCK EQ UATION 

Many equations describing evolution of physical, chemical, biological and social processes 
are written as mean-field equations for averages of quantities, which intrinsically are ran- 
dom processes to some extent. To take fluctuations around the averages into account, a 
straightforward way to proceed is to introduce a source of randomness directly in the mean- 
field equation. Then the quantities solved from the mean-field equations become stochastic 
processes depending on coordinate variables, i.e. stochastic fields. 

The paradigmatic example of this procedure is the Langevin equation for random walk, 
which describes the position r of a test particle subject to random force 77 

dr 

Tt =ri 

Here, the random force is of zero mean and uncorrelated in time (white noise), i.e. 

fai(*)»&(0> = i>M(*-0- 

Strictly speaking, this standard physical formulation is mathematically inconsistent, which 
gives rise to inevitable ambiguities in the case of multiplicative noise (when the noise term 
is multiplied by a function of the random position). 

Consider, for instance, the average distance between two points of the path of the random 
walk, i.e. (d is the dimension of space) - a hint to the property that the Brownian path, 
although continuous, is not different iable anywhere 

Increments of the Brownian path 



0)= / dtr](t) 

J to 



W(t) - W(t 

constitute an extremely important random process, the Wiener process, whose conditional 
probability density is the Gaussian 

V(W t\Wn tn) = - -(W-W ) 2 /2(t-to) 
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In particular, The Wiener process is the basis of the mathematically consistent definition of 
the Langevin equation (the stochastic differential equation, SDE). 

Critical dynamics. In the Landau theory of phase transitions the dynamics of the or- 
der parameter (p near equilibrium are described by the kinetic equation [time-dependent 
Ginzburg-Landau (TDGL) equation] 



;n -r^-W + ^+g^J- (63) 

In linear response theory, dynamics of fluctuations near equilibrium [28] are often described 
by kinetic equations similar to (63), but with a random noise term added to the right-hand 
side. 

In a more generic setup, standard models of critical dynamics are based on nonlinear 
Langevin equations 

where H is the effective equilibrium Hamiltonian. For the random source a suitable Gaus- 
sian distribution is assumed in which the correlation function is determined through the 
connection to the static equilibrium (fluctuation-dissipation theorem). For instance, model 
A for the non-conserved order parameter [29] is described by the SDE obtained from 63 by 
the addition of a white-noise field to the right-hand side. 

In reaction kinetics and population dynamics the simplest kinetic description of the dy- 
namics of the average particle numbers is given by the rate equation. The rate equation is 
a deterministic differential equation for average particle numbers in a homogeneous system, 
therefore it does not take into account boundary conditions, spatial inhomogeneities and 
randomness in the individual reaction events. Spatial dependence is often accounted for by 
a diffusion term, which gives rise to models of diffusion-limited reactions (DLR). 

As a simple example, consider the coagulation reaction A + A — > A. The diffusion-limited 
rate equation for the concentration if of the compound A is 

— = DV <p - kip , 

where k is the rate constant. 

The most straightforward way to take into account various effects of randomness is to 
add a random source and sink term to the rate equation: 

^ = DVV " V + / • (65) 
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This is a nonlinear Langevin equation for the field ip. Physically, in the case of concentration 
<P > o. 

There is an important difference between the reaction models and the critical dynamics: 
in the latter, deviations of the fluctuating order parameter from the (usually zero) mean 
may physically be of any sign (or direction). In particular, deviations from the equilibrium 
value are always allowed. In the reaction there is often an absorbing steady state, which 
does not permit fluctuations therefrom: once the system arrives at the absorbing state, it 
stays there forever. In particular, if the empty state is an absorbing state of the reaction, 
the the random source should be introduced multiplied by a factor vanishing in the limit 
ip — > to prevent the system returning from the absorbing state by the noise. The simplest 
choice yields 

^ = DV 2 p - kip 2 + ftp 

instead of (65). This is an equation with a multiplicative noise. 

1.4. Multiplicative noise. 

Consider the Langevin equation with the multiplicative noise of generic form 

^ = V{p) + fb(p) := -Kp + U{p) + fb(<p) , (66) 

where / is (usually) a Gaussian random field with zero mean and the white-in-time correla- 
tion function 

(f(t, x)f{t\ x')) = D(x - x>) = 5(t - t')D(x - x') , (67) 

where the shorthand notation x = (t,x) has been used. In (66), b(ip) is a functional of <p 
and U(ip) is a nonlinear functional of ip. Both functionals are time-local, i.e. depend only 
on the current times instant of the SDE. 

The Langevin equation with white-in-time noise / is mathematically inconsistent, because 
the time integral of the noise J fdt is a Wiener process which not differentiable anywhere as 
a function of time. 

This problem may be approached by starting with the set of correlation functions con- 
sisting of a 5 sequence in time, i.e. 

(/(*, x)f(t', x')) = D(t, x; t', x') ^ 8(t - t')D(x, x') (68) 

and passing to the white-noise limit at a later stage. >From the mathematical point of 
view, this treatment gives rise to the solution of the stochastic differential equation (66) 
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in the Stratonovich sense [30]. Physically, this is often the most natural way to approach 
the white-noise case. However, technically the Stratonovich interpretation gives rise to a 
rather complicated treatment and the SDE is most often used in the Ito interpretation in 
mathematical analyses. 

1.5. Fokker-Planck equation. 

Recall that the point of introducing of the SDE with white noise is to avoid dealing 
with the limit (68) explicitly. To this end, instead of using the mathematically problematic, 
although physically transparent, Langevin equation the stochastic problem (66), (67) may 
be equivalently stated in terms of the Fokker-Planck equation (FPE), which is an equation 
for both the conditional probability density p (ip, t\ipo, to) and the probability density p (ip, t) 
of the variable (p. Recall that both the master equation and the Fokker-Planck equation 
are special cases of the generic (forward) Kolmogorov equation and thus the two problems 
discussed here are closely related. The simple way to demonstrate this equivalence uses 
rules of Ito calculus [30], which is beyond the scope of the present treatment. Therefore, 
only the correspondence between the quantities specifying the stochastic problem in both 
approaches will be quoted here. The main advantage of the Fokker-Planck equation is that 
the equation itself is completely well-defined partial differential (or functional-differential for 
field variables) equation. The ambiguity of the Langevin problem shows in that the FPE is 
different for different interpretations of the SDE. 

For simplicity of notation, consider zero-dimensional field theory. The Fokker-Planck 
equation for the conditional probability density p (ip, t\ipo, to) m the case of the Ito equation 
is 

d d 

—p(ip,t\ip ,t ) = -—{[-K(p + U(ip)]p(p,t\p ,t )} 

+ \w ^ Db ^ p (v? ' to)] • (69) 

If the SDE (66) is interpreted in the Stratonovich sense, the FPE is 
d d 

—p(ip,t\ip ,t ) = -—{[-Kp + U(p)]p(p,t\p ,t )} 

+ ^{ %) ^ [jD%)p(v? ' t|v?0 ' to)] } • (70) 
The conditional probability density p(p,t\p ,t ) is the fundamental solution of the FPE 
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(69) or (70), i.e. 

p(ip,t \ip ,t ) = 5(ip - (f ) . 

Contractions are not quite obvious, when the random variable has several components. For 
instance, the Fokker-Planck equation in the Ito form becomes 

d d 

—p(ip,t\ip ,t ) = -— {[-Kijipj + Ui((p)]p((p,t\(po,t )} 

1 d 2 

+ 2 dip -dip - ^ hik ^ Dklh i l ^ p (fP't&O'to)] 

for a multi-component variable ip^. 

The Fokker-Planck equation may be regarded as the Schrodinger equation with imaginary 
time. Using this analogy, the solution of the FPE as well as calculation of expectation values 
may be represented in a way analogous to quantum field theory [31]. Construction with the 
FPE as the starting point gives rise to the famous Martin-Siggia-Rose solution of the SDE 
[20], but avoids ambiguities inherent in the SDE (they have been fixed by the choice of the 
FPE). 

Consider, for definiteness, the Fokker-Planck equation (69) corresponding to the Ito in- 
terpretation of the Langevin equation (66). Introduce - in analogy with Dirac's notation in 
quantum mechanics - the state vector | p t ) according to the following representation of the 
PDF 

p(<p,t) = (<p\pt), 

which is the solution of the FPE (69) with the initial condition p(tp, 0) = Po('p)- To construct 
the evolution operator for the state vector, introduce momentum and coordinate operators 
in the manner of quantum mechanics by relations 

d 

In these terms, the FPE for the PDF gives rise to the evolution equation for the state vector 
in the form 

J^Pt) = -H\pt), 

where the "Hamilton" operator for the FPE corresponding to the Ito interpretation of the 
SDE assumes, according to (69), the form 

H = -7T [-Kip + U{ip)\ - -n 2 b(0)Db(0) . (71) 

2 
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Note that, contrary to quantum mechanics, there is no ordering ambiguity in the construction 
of the Hamilton operator here. In this notation, the conditional PDF may be expressed as 
the matrix element 

p((p,t\(p ,to) = (if 

and the functional-integral representation may be constructed in the same fashion as in the 
master-equation case. 

2. FIELD-THEORETIC STUDY OF REACTION PROCESS A + A ->■ 

2. 1 . FIELD-THEORETIC MODEL OF ANNIHILATION PROCESS 

Let us study anomalous kinetics of the general type of the irreversible single-species 
annihilation reaction 

A + A 0, (73) 

with the unrenormalized (mean field) rate constant K . The first step of the Doi approach 
[8] (see also [23]) consists of the introduction of the creation and annihilation operators ^ 
and ip and the vacuum state |0) satisfying the usual bosonic commutation relations 

^(x),^(x')] = 5(x-x'), 
[^(x),^(x')] = [^t(x),^V)] = o, 

^(x)|0)=0,(0|^(x) = 0,(0|0) = l. (74) 

Let P({rii},t) be the joint probability density function (PDF) for observing particles at 
positions Xj. The information about the macroscopic state of the classical many-particle 
system may be transferred into the state vector | <&(£)) defined as the sum over all occupation 
numbers 

l^)> = E P (W'*)l{ni}>, (75) 
{«»} 

where the basis vectors are defined as 

IM>=n^( x *)P|0>. (76) 

i 

The whole set of coupled partial differential equations for the PDFs may be rewritten in the 
compact form of a master equation [7, 8] 

£|*(f)> = -H\*(t)h (77) 



-H(t-t ) 



(72) 
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where H = Ha + H D + H R and for the annihilation process A + A — > under consideration 

Ha = J rfx^V[v(x,t)^(x)], 
# D = -D J dx^V^x), 

# R = X D J dx(^) V, (78) 

corresponding to the advection, diffusion and reaction part [15]. Due to dimensional reasons 
we have extracted the diffusion constant D from the rate constant K = XoD . The mean 
of a physical quantity A(t) may be expressed [15] - with the use of the notation of Sec. 1 - 
as the vacuum expectation value 

(A(t)) = (0\T[A N ({l^(t)}) e~ ^ "' dt+n ^ d ^^' 0) ]\0}. (79) 

Here, the interaction operator is defined as H\ — H' — H' and the substitution ^ — >■ ift + 1 : 
H' = H(ipi + is understood. The field operators (74) have been replaced by the time- 
dependent operators of the interaction representation 

ip^t, x) = e H ° t il) i xeT H ° t , ip(t, x) = e^V(x)e~^o*. 

In this formulation - assuming the Poisson distribution as the initial condition - the average 
number density can be computed via the expression 

ra(t,x) = (Ol^x^V^^lO). (80) 

The expectation value of the time-ordered product in (79) can be cast [38] into the form of 
a functional integral over scalar fields V^( x j an d ^(x, t): 

(A(t)} = J V^A N {{l^{t)})e s \ (81) 
where the action S\ for the annihilation reaction A + A —> is 

dt J dx^dtip + ^V(v^) - D ^ f VV + 
AoA)^ 1 " + (V> f ) V + n y dx^(x, 0)}. (82) 

In order to analyze the effect of velocity fluctuations on the reaction process we average the 
expectation value (81) over the random velocity field v. The most realistic description of 
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the velocity field v(x) is based on the use of the stochastic Navier-Stokes equation. Due to 
the incompressibility conditions V.v = and V.f = imposed on the velocity field v and 
the random-force field f v it is possible to eliminate pressure from the Navier-Stokes equation 
and hence it is sufficient to consider only the transverse components 

d t v + P(v.V)v - z/ V 2 v = P. (83) 

Here, u is the molecular kinematic viscosity, -Pjj(k) = Sij — kikj/k 2 is the transverse pro- 
jection operator and k — |k| is the norm of the wave vector k. Here and below we use 
the subscript "0" for all "bare" parameters to distinguish them from their renormalized 
counterparts, which will appear during the renormalization procedure. 

The large-scale random force per unit mass f is assumed to be a Gaussian random variable 
with zero mean and the following correlation function 

(/m(Xi,*i)/„(x2,t 2 )) = <*(*1 ~*2) X 

/ J^PU^)df(k)^-^, (84) 

where the kernel function is chosen in the form 

df(k) = g 10 is 3 k*- d ~ 2t + g 20 u 3 k 2 . (85) 

The nonlocal term is often used to generate the turbulent velocity field with Kolmogorov's 
scaling [18, 19, 39]. This case is achieved by setting e = 2. The local term g 2 o^ok 2 has 
been added not only because of renormalization reasons but has also an important physical 
meaning. Such a term in the force correlation function describes generation of thermal fluc- 
tuations of the velocity field near equilibrium [17] and thus can mimic the usual environment 
in which chemical reactions take place. 

Averaging (61) over the random velocity field v is done with the "weight" functional 
W 2 = e S2 , where S 2 is the effective action for the advecting velocity field 

S 2 = iy dtdxdx' v(x,t).v(x',t)rf / (|x-x'|) + J dtdx v.[-<9 t v- (v.V)v + z/ V 2 v]. 

(86) 

With the use of the complete weight functional 



W = e s i+ s 2 



(87) 
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the expectation value of any desirable physical quantity is possible to be calculated. 

Actions (82) and (86) for the studied model are written in the form convenient for the 
use of the standard Feynman diagrammatic technique. There is, however, one delicate point 
which should be noted. In the pure reaction case we were dealing with the Cauchy problem 
and, consequently, all time integrals in the dynamic action were written over the positive 
time axis. When the random drift is included, it is rather natural and technically much 
simpler to regard the drift as a stationary random process given on the whole time axis. 
Then the integration over the fields ip and ift should include negative time arguments as 
well for consistency. In the pure reaction part this does not make any difference, because the 
retarded propagators render the time integrations at interaction vertices insensitive to the 
lower limit of integration. In the velocity part this leads to immense technical simplification, 
because the translation invariance with respect to time is preserved in construction of the 
perturbation expansion with the only exception of the initial conditions for the fields ip and 
ipi . In the case of Poisson initial distribution the initial condition is expressed as a linear 
term in the action and, in particular, does not affect renormalization of the model. 

Thus, henceforth we assume all time integrals in the dynamic action to be taken over the 
whole time axis and may use the standard Fourier representation to express the Feynman 
rules for the model. The inverse matrix of the quadratic part of the actions determines the 
form of the bare propagators. It is easily seen that the studied model contains three different 
types of propagators in Figure 1. In the momentum-frequency representation they are given 
as 

A^V,,k)= , ] (88) 
and in the momentum-time representation as 

(t, k) = 9(t) exp(-D k 2 t). (89) 

The vertex factor 

5 m V($) 

V m (x 1 ,x 2 ,...,x m ;$)= — — r (90) 

d^(x 1 )d^(x 2 ) . ..5®(x m ) 

is associated to each interaction vertex of a Feynman graph. Here, $ could be any member 

from the set of all fields {ip^, ip, v, v}. The interaction vertices from action (86) describe 

interactions between and it may be rewritten in a technically more convenient form 

- / dte v(v9)v = - / dto vvAv, = / 0«H»<, 
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where the incompressibility condition d(Vi = and partial integration method have been 
used. We assume that the velocity fields fall off rapidly for |x| — > oo. Rewriting this 
functional in the symmetric form ViVijiVjVi/2, it is easy to find the explicit form for the 
corresponding vertex factor in the momentum space 



Here, the momentum k is flowing into the vertex through the field v. The advecting term 
from the action (82) can be similarly presented as 



where the momentum k represents the momentum flowing into the vertex through the field 
ift . The vertices vvv and tp^ipv are depicted graphically in Figure 2. The two reaction 
vertices derived from the functional (82) according to the definition (90) are depicted in 
Figure 3 and physically describe the density fluctuations of the reactant particles. 

It should be stressed that in our model there is no influence of the reactants on the velocity 
field itself. Therefore, the model given by actions (82) and (86) may be characterized as a 
model for the advection of the "passive" chemically active admixture. 



The functional formulation provides a theoretical framework suitable for applying meth- 
ods of quantum field theory. Using RG methods it is possible to determine the IR asymptotic 
(large spatial and time scales) behaviour of the correlation functions. First of all a proper 
renormalization procedure is needed for the elimination of ultraviolet (UV) divergences. 
There are various renormalization prescriptions applicable for such a task, each with its own 
advantages. To most popular belong the Pauli-Villars, lattice and dimensional regulariza- 
tion [21]. In what follows we will employ the modified minimal subtraction (MS) scheme. 
Strictly speaking, in the analytic renormalization there is no consistent MS scheme. What 
we mean here, is the ray scheme [36], in which the two regularizing parameters e, A (e has 



Viji = i(kj6u + kiSij). 



(91) 





(92) 




2.2. UV RENORMALIZATION 
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been introduced in (85) and 2 A = d — 2 was introduced in (5)) are taken proportional to 
each other: A = £e, where the coefficient £ is arbitrary but fixed. In this case, only one 
independent small parameter, say, e remains and the notion of minimal subtraction becomes 
meaningful. UV divergences manifest themselves in the form of poles in the small expansion 
parameter and the minimal subtraction scheme is characterized by discarding all finite parts 
of the Feynman graphs in the calculation of the renormalization constants. In the modified 
scheme, as usual, certain geometric factors are not expanded in e, however. This is the 
content of the MS scheme used in our analysis. 

In order to apply the dimensional regularization for the evaluation of renormalization 
constants, an analysis of possible superficial divergences has to be performed. For the power 
counting in the actions (82) and (86) we use the scheme [18], in which to each quantity Q 
two canonical dimensions are assigned, one with respect to the wave number dq and the 
other to the frequency dq. The normalization for these dimensions is 

dZ = -d» = l, d k k = -d k x = l, d£ = d£ = 0. (94) 

The canonical (engineering) dimensions for fields and parameters of the model are derived 
from the condition for action to be a scale-invariant quantity, i.e. to have a zero canonical 
dimension. 

The quadratic part of the action (82) determines only the canonical dimension of the 
quadratic product ip^ip. In order to keep both terms in the nonlinear part of the action 

A A) j dtdx[2^ + (V> f ) V, (95) 

the field tp^ must be dimensionless. If the field ift has a positive canonical dimension, 
which is the case for d > 2, then the quartic term should be discarded as irrelevant by the 
power counting. The action with the cubic term only, however, does not generate any loop 
integrals corresponding to the density fluctuations and thus is uninteresting for the analysis 
of fluctuation effects in the space dimension d = 2. 

Using the normalization choice (94) we are able to obtain the canonical dimensions for all 
the fields and parameters in the rf-dimensional space. The results are summarized in Table 
1. 

Here, dq = dq + 2dq is the total canonical dimension and it is determined from the 
condition that the parabolic differential operator of the diffusion and Navier-Stokes equation 
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scale uniformly under the simultaneous momentum and frequency dilatation k — > fik, u — > 
fi 2 u. 

The model is logarithmic when all coupling constants g w , g 2 o, \o vanish simultaneously. From 
Table 1 it follows that this situation occurs for the choice e = A = 0. The UV divergences 
have the form of poles in various linear combinations of e and A. The total canonical 
dimension of an arbitrary one-particle irreducible Green (1PI) function r = (<&.. . $)i_i r is 
given by the relation d T = d + 2 — N^d®, where Nq> = {N^t,N^, N v , N^} are the numbers 
of corresponding external fields. The statistical averaging (. . .) means averaging over all 
possible realizations of fields v, v, ip satisfying appropriate boundary conditions with the 
use of the complete weight functional (87). Superficial UV divergences may be present only in 
those T functions for which dr is a non-negative integer. The superficial degree of divergence 
for a 1PI Green function T is 

d T = 4 - N v - Nt - 2AV (96) 

However, the real degree of divergence 5 T is smaller, because of the structure of the interaction 
vertex (91), which allows for factoring out the differential operator d to each external line v. 
The real divergence exponent 5r may then be expressed as 

5 r = d r - N s = 4 - N v - 2N^ - 2N i> (97) 

Although the canonical dimension for the field ip^ is zero, there is no proliferation of superfi- 
cial divergent graphs with arbitrary number of external ijy legs. This is due to the fact that 
n^t < rty, which may be established by a straightforward analysis of the graphs [7]. Brief 
analysis shows that the UV divergences are expected only for the 1PI Green functions listed 
in the Table 2. 

This theoretical analysis leads to the following renormalization of parameters g , D and 

u 

9i = 9wH~ 2ez l, 92 = g 2 oV 2A ZfZ- 1 , u = u^Z x Z~ x , 

A = \ /j 2A Z 2 Z^ 1 , u = u Z{\ D = D Z 2 ~ 1 , (98) 

where /j, is the reference mass scale in the MS scheme [21] and we have introduced the 
inverse Prandtl number u — D/vior convenience. It represents the ratio between diffusion 
and viscosity forces in a liquid. In terms of introduced renormalized parameters the total 
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renormalized action for the annihilation reaction in a fluctuating velocity field is 

S R = J dxdt^d^ + ^V(vV') - ui/Z 2 vV + \uufj,- 2A Z 4 [2^ + (V> f ) V + 
n | rfx^t( x , 0) - ^vfoi/V^-V 2 ) 1 -^ " <^V 2A ^V 2 ]v + 
v.[^v + (v.V)v - z/ZiV 2 v]|. 



(99) 



The renormalization constants Zj,i = 1,2,3,4 are to be calculated perturbatively through 
the calculation of the UV divergent parts of the 1PI functions T^^, T^^, r^tp^, and 
Interaction terms corresponding to these functions have to be added to the original 
action S — S± + S 2 with the aim to ensure UV finiteness of all Green functions generated by 
the renormalized action Sr. At this stage the main goal is to calculate the renormalization 
constants Z i: % — 1,2, 3, 4. 

The singularities in various Green functions will be realized in the form of poles in e and A 
and their linear combinations such as 2e + A or e — A. Recall that for the consistency of the 
MS scheme it is necessary that the ratio 

i = ^ (100) 

is a finite real number. It should be noted that the graphs corresponding to T^^2 and 
r^t)2^2 differ only by one external vertex and thus give rise to equal renormalization of the 
rate constant AoA)- Therefore, in what follows, we will always consider the function T^t^a. 
In order to calculate the renormalization constants Z 2 and Z 4 we proceed according to the 
general scheme suggested in [36]. We require the fulfillment of UV finiteness (i.e. finite limit 
when e, A — > ) of the 1PI functions r^ty,|u>=o an d r^t^a | w =o- Because the divergent part 
of the Feynman graphs should not depend on the value of uj, we have adopted the simplest 
choice oj = 0. It is convenient to introduce the dimensionless expansion variables of the 
perturbation theory as 

_ dioSd _ 92oSd _ ^oSd /iai\ 

"10 = -^, "20 = ^a, «30 = ^a, (101) 

where Sd is the surface area of the unit sphere in d— dimensional space, p is the total mo- 
mentum flowing into the Feynman diagram and Sd = Sd/(2n) d . For brevity, in the following 
we use the abbreviation g = g Q Sd for the parameters {#io, #20, Ao} or their renormalized 
counterparts, respectively. Next we demonstrate the perturbation series for the 1PI Green 
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functions to the second order approximation. The perturbative expansion for may be 
written as 



tpt^i \ u)=0 



D oP 2 



ni+n2=2 



1+ E «io«20 2 7^ n2) (rf,«o) 



ni+n 2 >l 



(102) 



where 7^,t^ are dimensionless coefficients which contain poles in e and A. Explicit dependence 
on the space dimension d and inverse Prandtl number u is emphasized. It is important to 
note that there are no terms in this sequence proportional to the expansion parameter a 30 . 
In terms of the renormalized parameters perturbative expansion for the Green function is 
(102) 



r 



Dp 2 



ni+n 2 =2 



ni,"2>0, 
?ii+n2>l 



(103) 



with the renormalized parameters a± = g~[s 2e Z^ and a 2 = 'g 2 s~ 2A Z 3 Z^ 3 in accordance with 
the relations (98) and (101), where s = fi/p. Here we would like to stress, that in order to 
get the correct expansion in e and A, one has to make replacement 



d ->■ 2 + 2A, m ->■ Z~ 1 Z 2 u, 



(104) 



in the arguments of 7^™ 2 ' ) - In the same way, the perturbation expansion series for the 
Green function Y^^2 is 

n.iH-n2+«3=2 



L=o — — 



1+ E "M^^^fatio) 

ni,ri2,n:j>0, 
ni+ri2+n 3 >l 



(105) 



where T^t^a are dimensionless coefficients resulting from calculation of Feynman graphs. 
Again by replacing the bare parameters with the renormalized counterparts the following 
series is obtained 

ni+ri 2 +«3=2 



AXDjj, 



-2A 



E 



(d,u) 



ni,n2,n3>0, 



(106) 



where the dimensionless parameter a 3 = Xs~ 2A Z 2 ~ 1 Z4 is introduced and the change (104) is 
understood. The perturbation series for the Green function T^tp^ has the same form, so 
we do not present it. 

Denoting by Z^ the contribution of the order g n , g = {gi, g 2 , A} the first order of renor- 
malization constants Z 2 and Z 4 may be calculated via equations 



■ 2 ~ L l9is +g 2 s y H 



(107) 
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r(l) 



r[?7-^ (1 ' ' 0) J_ TJ-q- 2 ^ ' 1,0 ) _L T„-2A (0,0,1)1 



(108) 



where C stands for the operation of extraction of the UV-divergent part (poles in e and A or 
their linear combination). In the MS scheme finite terms are discarded, so we do not need 
to take care of them. At the second order the term for Z 2 can be schematically written as 

u 



,(2) 



= C 



9is 2t 
l + u 



uZ { 2 l) + {u + 2)Z\ 



(1) \ (1,0) 2A 

h,,u -92S 



• z? + ^z? 

l + u 2 l+u 1 



7 (1) \ (0,1) 2 4, (2,0) Js^A^lU) 2 -4A f0,2) 

z 3 1 7^ +9i s y H + gig 2 s +92 s y H 



(109) 



The two-loop graphs that contribute to the calculation of Z 2 are represented by the graphs 
depicted in Figure 4. For the renormalization constants Z 4 we have the expression 



7 (2) _ n y 2( e -A) fl,0,l) y -4A (0,1,1) ,T 2 -4A (0,0,2) 2e J 1,0,0)/ oyWu 

z 4 - -H9i* s 7(^+)2^,2 + 92*s 7^ + A s 7^2 + g x s 7^ J + 



fcs-^'W " 3Z«) + ^"-7^(2^ - Z«)]. 



(110) 



The two-loop graphs that contribute to the calculation of Z4 are represented by the graphs 
depicted in Figure 5. >From these expressions the renormalization constants Z 2 and Z 4 can 
be calculated in the form 



9i 



+ 



92 



8u(l + u)e 8w(l + w)A 
BiigT 2 ^ B 2 292 2 _ Bi 2 gig 2 
~^ + ^A~ + 7^A~' 



+ ^ll^l 2 + ^22#2 2 + Al29l92 + 



A 2 



eA 



111) 



1 - 



A 



9i A 



92^ , A 2 



2A 16u(l + u) (e - A)A 32w(l + w)A 2 4A 2 
5-iA g 2 X\ , , 



(112) 



e-A A, 

The lengthy expressions for the coefficient functions Aij(£,u),Bij(£,u) and C(u, £) can be 
found in appendix A. 

In a similar way we obtain renormalization constants Z\ and Z 3 [36] from condition of 
the UV finiteness for the 1PI Green functions |^=o and r ss | w=0 . The perturbation series 
for Tvv can be written as 

ni+ri2=2 



vv\ui=0 



_1 4. ^ n 2 Ani,n 2 ) 

1 ^ 2-^i 10 20 



(d) 



ni,ri2>0, 
ni+n2>l 



(113) 
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and for as 



VV \U!=0 



3 2-2A-2e 



ni+U2=2 



ni>0,ri2>— 1, 
ni+n2>l 



(114) 



>From the definition of the projection operator Pf- it is easy to see, that after contracting 
indices i and j we are left with the constant d—1. Hence, rewriting perturbations series for 
F iv and r cc in the renormalized variables (98) and contracting indices i and j we get 

r. 



vv \u)=0 



1.1+712=2 



up 2 (d- l) 



ni,7i2>0, 
ni+n 2 >l 



(115) 



|u=0 



"1+12=2 

_ 9l 2 £ +2A , y , 7 n. ni n. ri2 o/( ni '' 12 ) 

(rf^W " S s + z * + z ' x „ ° 2 7 « 

ni+n2>l 
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By the same algorithm as described above in detail for the calculation Zi and Z4 explicit 
expressions for the renormalization constants Z\ and Z3 are obtained. The results for them 
in the MS scheme can be found in [36]. 

2.3. IR STABLE FIXED POINTS AND SCALING REGIMES 



The coefficient functions of the RG differential operator for the Green functions 



Drg = A* 



d_ 

dp, 



d_ 

dp, 



9i 



_d_ 

'dgi 



d_ 

dv : 



where the bare parameters are denoted with the subscript "0", are defined as 



71 = pt- 



d In Z 1 



dpi 



A = p- 



dgi 
dp, 



(117) 



(118) 



with the charges gi = {gi, g2, u, A}. From this definition and the renormalization relations 
(98) it follows that 



(3 gi = 9l (-2e + 3 7 i) , (3 g2 = g 2 (2A + 3 7l - 73), 
/3 A = A(2A- 7 4 + 7 2 ), /9« = «(7i-72), 

where the anomalous dimensions 7^ (a = 2,3,4) are defined as 

<91nZ n 



7a = H- 



dp, 



(119) 



(120) 
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We are interested in the IR asymptotics of small momentum p and frequencies oj of the renor- 
malized functions or, equivalently, large relative distances and time differences in the (t, x) 
representation. Such a behaviour is governed by the IR-stable fixed point g* = (#*, g%,u* , A*), 
which are determined as zeroes of the /3 functions /3(g*) = 0. The fixed point g* is IR stable, 
if real parts of all eigenvalues of the matrix Uij = d(3i/dgj\ g=g * are strictly positive. From the 
knowledge of renormalization constants Z 2 and Z 4 (111), (112) and definitions (118), (120) it 
is possible to calculate anomalous dimensions 72 and 74 

^ = Tirrh ~ 4B ^ 2 + 4£?22 ^ 2 " 2S i2^> (121) 

Au(l + u) 

74 = _A + A(^r + ^)C(«,0- (122) 

A straightforward calculation shows that higher order poles cancel each other, so that the 
anomalous dimensions 72 and 74 are finite. For completeness we quote also anomalous 
dimensions 71 and 73 [36] to the same order 

gi + gi , (4£ + 3) _ 2 5^ + 3 R ( _^_, 2 , 10 ,, 

71 = -l6~ + ^TO 91 + " 256^ + ^ ' (123) 

= (gT + gl) 2 _ e(13 + 19Q ^I 3 34e + 19 + 6^ _ 2 _ 3^ 13 + 31^ 

73 16^ 1024(2 + f) 02 512(2 + ^ 512 1024 ^ 2 

i^(E±M! (124) 

256 ^ 

where the value i? = —0.168 is a result from numerical integration. Zeroes of the beta 
functions (119) determine possible IR behaviour of the model. There are four IR stable fixed 
points and one IR unstable fixed point. In this section we present them with their regions 
of stability. 

(i) The trivial (Gaussian) fixed point 

g? = g?=T=0, (125) 

with no restrictions on the inverse Prandtl number u. The Gaussian fixed point is stable, 
when 

e<0, A>0. (126) 

and physically corresponds to the case, when the mean-field solution is valid and fluctuation 
effects negligible. 
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(ii) The short-range (thermal) fixed point 

m* = 0, gl* = -16A + 8(l + 2i?)A 2 , 

= VTf- 1 _ L12146A 
2 

A* = -A + 2.64375), (127) 

at which local correlations of the random force dominate over the long-range correlations. 
This fixed point has the following basin of attraction 

a 2i2 - 1 A 9 3 A 2 . . 

A — A < 0, 2e + 3A - — - < 0, (128) 

A + ^A 2 < 0, A + 0.4529Ae < (129) 
2 

and corresponds to anomalous decay faster than that due to density fluctuations only but 
slower than the mean-field decay. 

(ra) The kinetic [37] fixed point with finite rate coefficient 

_ 32 e (2e + 3A) 

92 = ~9~ ATe ' 

, 7T7-1 
M = g + M i(0 e > 

A* = -|(c + 3A) + i-(3A + e)(Qe - A), (130) 
Here Q = 1.64375. The fixed point (130) is stable, when inequalities 

^±>0, e>0, -jje<A<-ie, (131) 

are fulfilled, where 

3 3 9 y 

4e(e + 3A)R - 6e 2 - 12eA - 9AM 



V9A 2 - 12Ae - 8e 2 

The decay rate controlled by this fixed point of the average number density is faster than the 
decay rate induced by dominant local force correlations, but still slower than the mean-field 
decay rate. 
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(iv) The kinetic fixed point with vanishing rate coefficient: 

32 e(2e + 3A) 2 
9i = T e + A +9A&, 

^ 32 6 . . . n 

«* = ^|zl + ^) c> T = 0. (133) 

This fixed point is stable, when the long-range correlations of the random force are dominant 

ft±>0, e>0, A>-^e, (134) 

o 

and corresponds to reaction kinetics with the normal (mean-field like) decay rate. 
(v) Driftless fixed point given by 

gl* = gi* = 0, u* not fixed, X = -2A, (135) 

with the following eigenvalues 

fii = -2e, tt 2 = -fl 4 = 2A, tt 3 = 0. (136) 

An analysis of the structure of the fixed points and the basins of attraction leads to 
the following physical picture of the effect of the random stirring on the reaction kinetics. 
Anomalous behaviour always emerges below two dimensions, when the local correlations are 
dominant in the spectrum of the random forcing [the short-range fixed point (ii)]. However, 
the random stirring gives rise to an effective reaction rate faster than the density-fluctuation 
induced reaction rate even in this case. The anomaly is present (but with still faster decay, 
see the next Section) also, when the long-range part of the forcing spectrum is effective, but 
the powerlike falloff of the correlations is fast [this regime is governed by the kinetic fixed 
point (Hi)]. Note that this is different from the case in which the divergenceless random 
velocity field is time-independent, in which case there is no fixed point with A* 7^ 0[12]. At 
slower spatial falloff of correlations, however, the anomalous reaction kinetics is replaced by 
a mean-field-like behaviour [this corresponds to the kinetic fixed point (iv)]. In particular, 
in dimensions d > 1 this is the situation for the value e = 2 which corresponds to the 
Kolmogorov spectrum of the velocity field in fully developed turbulence. Thus, long-range 
correlated forcing gives rise to a random velocity field, which tends to suppress the effect of 
density fluctuations on the reaction kinetics below two dimensions. For better illustration, 
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regions of stability for fixed points (i)-(iv) are depicted in Fig. 6. Wee see that in contrast to 
the one-loop approximation [15], overlap (dashed region) between regions of stability of fixed 
points (ii) and (Hi) is observed. It is a common situation in the perturbative RG approach 
that higher order terms lead to either gap or overlap between neighbouring stability regions. 
The physical realization of the large-scale behaviour then depends on the initial state of the 
system. 

2.4. LONG-TIME ASYMPTOTICS OF NUMBER DENSITY 

Since the renormalization and calculation of the fixed points of the RG are carried out 
at two-loop level, we are able to find the first two terms of the e, A expansion of the 
average number density, which corresponds to solving the stationarity equations at the one- 
loop level. The simplest way to find the average number density is to calculate it from 
the stationarity condition of the functional Legendre transform [38] (which is often called 
the effective action) of the generating functional obtained by replacing the unrenormalized 
action by the renormalized one in the weight functional. This is a convenient way to avoid 
any summing procedures used [7] to take into account the higher-order terms in the initial 
number density n . We are interested in the solution for the number density, therefore we put 
the expectation values of the fields v and v equal to zero at the outset (but retain, of course, 
the propagator and the correlation function). Therefore, at the second-order approximation 
the effective renormalized action for this model is 





X 

(137) 

where S\ is the action (82) (within our convention Si = in the effective action) and graphs 
are shown together with their symmetry coefficients. The slashed wavy line corresponds 
to the field ip^ and the single wavy line to the field ifi. The stationarity equations for the 
variational functional 
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give rise to the equations 

d t i\) = uuZ 2 V 2 ip - 2Xuufi- 2A Z 4 (1 + ^) i? 

oo 

+4wVA 2 /i~ 4A J dt' j dy (A^) 2 (t -t',x- y)^ 2 (t', y) 



oo 


oo 

-jdt'j dy A- (t -t',x- y) ^A^ f (t -t',x- y)^(t>, y) + . . . , (139) 



oo 

d_ 

dx 



o 



-d t ^ = uuZ 2 V V - 2\uufi~ 2A Z 4 2^+(^) V 

oo 

+8«VaV 4A J dt' j dy(A^) 2 (t' -t,y- x)^(t', y)^(t, x) 
o 

oo 

+4 M VAV 4A j "dt! j 'dy(AW 'fit' -t,y - x) x [^(t',y)} 2 i/i(t,x) + 
o 

oo 

J dt' I dyA™(t' -t,y- x)^-A^\t' -t,y- </) + ••■ (140) 

o J 
In (139) and (140), in the integral terms it is sufficient to put all renormalization con- 
stants equal to unity. Substituting the solution ^ — of (140) into (139) we arrive at the 
fluctuation-amended rate equation in the form 

d t ij = uuZ 2 V 2 ^ - 2\uufi~ 2A Z 4 ij 2 (141) 

oo 

+4 M 2 z/V 4A A 2 j dt' j dy {A^) 2 {t -t',x- y)ij 2 (t', y) 
o 

oo 

+ ^Jdt' jdy A»J (t -t',x- y) {t _ ^ x _ „) + ..., (142) 

o J 
This is a nonlinear partial integro-differential equation, whose explicit solution is not known. 

It is readily seen that for a homogeneous solution the term resulting from the third graph in 

(137) vanishes and hence the influence of the velocity field on the homogeneous annihilation 

process would be only through the renormalization of the coefficients A and D. However, in 

case of a nonuniform density field ip the effect of velocity fluctuations is explicit in (141). 

Such a solution can be most probably found only numerically. 



37 



To arrive at an analytic solution, we restrict ourselves to the homogeneous density n(t) = 
(ip(t))i which can be identified with the expression (62). In this case the last term in (141) 
vanishes together with the Laplace operator term and the remaining coordinate integral may 
be calculated explicitly . The propagator is the diffusion kernel of the renormalized model 
(we consider first the system in the general space dimension d) 



A^ T (t - t', x) = ^ '—^ exp 

[Anuu(t - t')f 2 



x 



4uv(t - f ) 



(143) 



As noted above, for calculation of the one-loop contribution it is sufficient to put the renor- 
malization constant Z2 = 1 in the propagator . Therefore, evaluation of the Gaussian 
coordinate integral in (141) yields 

[8nuu{t — 1')\ 1 

and we arrive at the ordinary integro-differential equation 



t 



^p- = -2\uv^Z,n\t) + A\ 2 u 2 u 2 fi- 4A [ dt' U ^ . (145) 

dt J [87TMZ/(t — t')] 1 



Spatial fluctuations in the particle density show in the integral term and affect rather heavily 
even the homogeneous solution. In particular, the integral in (145) diverges at the upper 
limit in space dimensions d > 2. This is a consequence of the UV divergences in the model 
above the critical dimension d c = 2 and near the critical dimension is remedied by the UV 
renormalization of the model. To see this, subtract and add the term n 2 (t) in the integrand 
to obtain 

t 

^P- = -2\uvfi- 2A Z,n 2 (t)+4\ 2 u 2 v 2 fi- 4A n 2 (t) [ — — 

t 

+4A WaT 4A / dt' ^)-n 2 (t) (14g) 
J [8nuu(t - t')f 2 

The last integral here is now convergent at least near two dimensions, provided the solution 
n(t) is a continuous function. This is definitely the case for the iterative solution constructed 
below. The divergence in the first integral in (146) may be explicitly calculated below two 
dimensions and is canceled - in the leading order in the parameter A = (d — 2)/2 - by the 
one- loop term of the renormalization constant Z4 (112). Expanding the right-hand side of 
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(146) in the parameter A = (d — 2)/2 to the next-to-leading order we arrive at the equation 
dn(t) 



dt 



= -2\uufi- 2A n 2 (t) + 2\uufi- 2A n 2 (t) j A [7 + In (2uufi 2 t)] j 

' 2 " 2A /^" 2(t ;';; 2( ^ (147) 



A 2 w/i~ 2A f^n 2 (t')-n 2 (t) 




without divergences near two dimensions. Here, the factor fi~ 2A has been retained intact in 
order not to spoil the consistency of scaling dimensions in different terms of the equation. In 
(147), 7 = 0.57721 is Euler's constant and we have considered the coupling constant A and 
the parameter A = (d — 2)/2 to be small parameters of the same order taking into account 
the magnitudes of the parameters in the basins of attraction of the fixed points of the RG. 
The leading-order approximation for n(t) is given by the first term on the right-hand side 
of (147) and it is readily seen that after substitution of this expression the integral term in 
(147) is of the order of A 3 and thus negligible in the present next-to-leading-order calculation. 
In this approximation, Equation (147) yields 

n(t) = , (MS) 

1 + 2Xuut 1 1 + [1 - 7 - In (2uu(jlH)] j fi- 2A n 

where n is the initial number density. 

Since the fields $ = {v , v, ip, ip^} are not renormalized, the renormalized connected Green 
functions Wr differ from the unrenormalized W = ($...$) [19] only by the choice of 
parameters and thus one may write 

W R (g,v,»,...) = W(g ,v ,...), (149) 

where g = {gw, #20, Uo, Ao} is the full set of the bare parameters and dots denotes all 
variables unaffected by the renormalization procedure. The independence of renormalization 
mass parameter /i is expressed by the equation fid^Wu = 0. Using this equation the RG 
equation for the mean particle number n(t) is readily obtained: 

(d d d \ 

»djl + Pa-Qg - 7lZ/ ^J n ^ U > n °' 9) = °- ( 15 °) 

We are interested in long-time behaviour of the system (t — > 00), therefore we trade the 

renormalization mass for the time variable. Canonical scale invariance yields relations [34] 

d d d \ 
fx— - 2z/— + dn — d J n(t, //, v, n , g) = 0, (151) 
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d d \ 

-t-^; + v-r^jn(t,n,v,n o ,g)=0, (152) 

where the first equation expresses scale invariance with respect to wave number and the 
second equation with respect to time. Eliminating partial derivatives with respect to the 
renormalization mass u and viscosity v we obtain the Callan-Symanzik equation for the 
mean particle number 



dn 



n(t,n, v,n ,g) = 



(153) 



To separate information given by the RG, consider the dimensionless normalized mean par- 
ticle number 



n 



n 



n 



— = $ [ufi t,\u—;,g . 



(154) 



For the asymptotic analysis, it is convenient to express the particle density in the combination 
used here. Solution of (153) by the method of characteristics yields 



n 



$ ( z//Tt, \u—,g) = $ ( z/u 2 r, \u 



no 



9 



(155) 



where r is the time scale. In Equation (155), g and Uq are the first integrals of the system 
of differential equations 



d__ 
l Jt 9 ~ 2- 7 i(S) 



d 



n 



— , t^-nn = d 

dl 2-7i( 



(156) 



Here g = {g 1 , g 2 , u, A} with initial conditions g\t= T = g and no|t= T = ^o- In particular, 



\un = \un Q 



exp 



74 (is 



(2-7i)s 



(157) 



The asymptotic expression of the integral on the right-hand side of (157) in the vicinity of 
the IR-stable fixed point g* is of the form 



74<is 



74 



In - + 



(74 -lt)ds 7| 



In - + c 4 (r) , (158) 



(2 - 7i)s t^oo 2 - 7* Vry 2 - 7* J (2 - 7l )s 2 - 7* V r 

T 

corrections to which vanish in the limit t — >■ 00. In (158) and henceforth, the notation 
7* — 7i (9*) nas been used. >From the point of view of the long-time asymptotic behaviour 
the next-to-leading term in (158) is an inessential constant. In the vicinity of the fixed point 
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where a shorthand notation y has been introduced for the long-time scaling of the normalized 
number density as well as the dimensional normalization constant 

and the decay exponent 

a = 1 + (160) 
1 ~ 7i 

The asymptotic behaviour of the normalized particle density is described by the scaling 
function f(x, y) 

$ (Vt, \u ^) ~ $ (Vr, C n y, g*) = / (Vr, C n y) . (161) 

The free parameters in the variables of the scaling function f(x,y) correspond to the choice 
of units of these variables, whereas the objective information is contained in the form of the 
scaling function [19, 34]. Here, it is convenient to use the explicit solution (148) to obtain 
the e, A expansion for the inverse h(x,y) = 1/ f(x,y) of the scaling function. We obtain the 
generic expression 

h(x, y) = j^—j = l + 2x 2 /jl + A_[l- 7 -ln (2u*x)} j , (162) 

the substitution in which of the various fixed-point values A* (at the leading order A* ~ 27rA*) 
and u* in the leading approximation yields the corresponding e, A expansions. 

Below, we list the scaling functions h(x, y) and the dynamic exponents a at the stable 
fixed points in the next-to-leading-order approximation. 

(i) At the trivial (Gaussian) fixed point (125) the mean-field behaviour takes place with 

h(x, y) = 1 + 2xy , 

a = l. (163) 

(ii) The thermal (short-range) fixed point (127) leads to scaling function and decay exponent 

h(x,y) = l + 2xy^l-^ [l - 7 - In (>/l7 - l) x] j , 

A A 2 . . 

« = l + y + T - (164) 

Here, the last coefficient is actually a result of numerical calculation, which in the standard 
accuracy of Mathematica is equal to 0.5. We have not been able to sort out this result ana- 
lytically, but think that most probably the coefficient of the A 2 term in the decay exponent 
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a in (164) really is \. 

(in) The kinetic fixed point with an anomalous reaction rate (130) corresponds to 

h(x,y) = l + 2xy!^l-^^ [l - 7 - In (y/vf - l) x] |, 

3A + e , . 

« = l + (165) 

with an exact value of the decay exponent. 

(iv) At the kinetic fixed point with mean- field-like reaction rate (133) we obtain 

h(x,y) = l + 2xy, 

a = l. (166) 

In the actual asymptotic expression corresponding to (161) the argument y — > C n y is different 
from that of the Gaussian fixed point. 

To complete the picture, we recapitulate - with a little bit more detail - the asymptotic 
behaviour of the number density in the physical space dimension d = 2 predicted within 
the present approach [15] (it turns out that for these conclusions the one-loop calculation is 
sufficient). On the ray e < 0, A = logarithmic corrections to the mean-field decay take 
place. The integral determining the asymptotic behaviour of the variable (157) yields in this 
case 

^ f'Ua.w, (167) 



(2 — 7i)s t^-oo 2 \t 
with corrections vanishing in the limit t — > oo. Therefore, in the vicinity of the fixed point 

3 "^~ A v(7) ln " /2 (?) <5 "" ?7e '" (168) 

The scaling function h is of the simple form 

h(x, y) — 1 + 2xy 

and gives rise to asymptotic decay slower than in the mean-field case by a logarithmic factor: 

In 1 / 2 (t/r) 

n ~ ~ — . 

2u\uC n t 

It is worth noting that this logarithmic slowing down is weaker than that brought about the 
density fluctuations only [23] and this change is produced even by the ubiquitous thermal 
fluctuations of the fluid, when the reaction is taking place in gaseous or liquid media. 
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On the open ray e > 0, A = the kinetic fixed point with mean-field-like reaction rate 
(133) is stable and the asymptotic behaviour is given by (166) regardless of the value of the 
falloff exponent of the random forcing in the Navier-Stokes equation. In particular, only the 
amplitude factor in the asymptotic decay rate in two dimensions is affected by the developed 
turbulent flow with Kolmogorov scaling, which corresponds to the value e — 2. This is in 
accord with the results obtained in the case of quenched solenoidal flow with long-range 
correlations [12, 13] as well as with the usual picture of having the maximal reaction rate in 
a well-mixed system. 

3. ROLE OF RANDOM SOURCES AND SINKS ON REACTION PROCESSES 

3.1. MASTER EQ UATION FOR RANDOM SO URCES AND SINKS 

We will consider the annihilation reaction A + A — > in a random drift field in a more 
general setup than in the previous parts. For this purpose we introduce random sources and 
sinks of the reacting particles in order to maintain a steady state in the system. In most 
cases this is carried out by including an additive noise term in the Langevin equation of 
the stochastic process as was done e.g. in (83) to have steady turbulent state. Since our 
analysis is based on the master equation, this is not quite appropriate here. Unfortunately, 
there is no unique way to introduce random sources in the master equation corresponding 
to the random noise of the mean-field (Langevin) description. We use the simplest choice, 
described in detail in [4] , which is equivalent to adding processes A — > X and Y — > A .to the 
whole reaction scheme. Here X and Y stand for particle baths of the sink and the source 
respectively. In a homogeneous system these reactions leads to the master equation 

dP ^' = n + V [P(t, n - 1) - P(t, nj\ + [{n + l)P(t, n + 1) - nP(t, nj\ . . . (169) 

where P{t, n) is the probability to find n particles at the time instant t in the system. The 
ellipsis in (169) represents terms describing the annihilation reaction, diffusion and advection 
in the system. In (169) fi + and yU_ are the reaction constants of the creation and annihilation 
reactions, respectively. The transition rate has been chosen proportional to the particle 
number n, which can be understood as consequence of independent processes A — > X and 
this choice also preserves the empty state as an absorbing state. In the transition rate for 
creation process V is the volume of the (for the time being) homogeneous system and will 
be important in passing to the continuum limit of the inhomogeneous system. The master 
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equation (169) gives rise to the reaction-rate equation 

^ = pL+ V-pL-(n) + ... (170) 

where (n) is the mean particle number. 

We recall that the basic idea of the Doi approach [8] is to rewrite the set of master 
equations for probability distributions of a stochastic problem in the form of a single kinetic 
equation for a state vector incorporating all probabilistic information about the system 
constructed in a suitable Fock space. The kinetic equation is defined by the Liouville operator 
acting in the Fock space and generated by the set of master equations. Although the basic 
procedure has been thoroughly exposed in the literature, the introduction of random sources 
and sinks of particles in the master equation has specific features, which should be presented 
in detail. Therefore, let us briefly recall the basic quantities and relations of the Doi approach. 
For simplicity, consider probabilities P(t,n) to find n particles at the time instant i on a 
fixed lattice site. Then the spatial dependence may be described by labeling the particle 
number by the coordinates of the lattice and introducing necessary sums and products over 
the lattice sites. The construction of corresponding Fock space was presented in Section. 1, 
namely equations (7-14). The set of master equations for a birth-death process may also be 
cast in the form of a single evolution equation for the state vector (14) without any explicit 
dependence on the occupation number 

^ = -tf(a + ,a)|$>. (171) 

Master equations (169) give rise to the following terms in the Hamilton operator 

H g (a + , a) = -fi + V (a+ -/)-//_(/- 6+) a , (172) 

where / is the identity operator. The expectation value of any function A(n) of the random 
particle number 

oo 

(A(t)) = J2A(n)P(t,n), (173) 

n=0 

may be expressed in the form of the functional integral over the functions a(t) and a{t) 

(A(t)) = J VaDa A N {1, a{t))e Sl , (174) 
where An(cl, a) is the normal form [38] of the operator A{a + a) and S\ is the dynamic action 

oo 

Si(a,a) = Jdt [-a(t)d t a(t) + fx + Va(t) - fi_a(t)a(t)] . . . (175) 
o 
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Only the generic time-derivative term and terms brought about by the random source model 
are expressed here explicitly, while the ellipsis stands for terms corresponding to other reac- 
tions and initial conditions. 

Let the transition rates /a± be the random functions uncorrelated in time with a probability 
distribution given in terms of the moments ([/,±) = E± >n . To keep the problem translation- 
invariant in time we assume stationary stochastic processes determined on the whole time 
axis. Therefore, henceforth all time integrals in the action functional and, correspondingly, 
in the perturbation expansion shall be taken over the whole time axis. At this point we 
also generalize the treatment to the case of a spatially inhomogeneous system and introduce 
a lattice subscript as the spatial argument, i.e. a(t) — > aj(t). In this case the volume V 
becomes the volume element attached to the lattice site. For simplicity, we replace the time 
integral with the integral sum dt — >• ^ Q At and assume that the transition rates at each 
time instant and lattice site n± !0t ,i are independent random variables. Then the average of the 
expectation value (174) over the distribution of random sources reduces to the calculation 
of the expectation value 

n( ex P ( ,a,iVa ati At - n_ ^ia a)i a ati At J \ . (176) 

a,i ' ' 

For each particular time instant and lattice (we assume that the moments of /i± are the same 
for all a and i and omit labels for brevity) this gives rise to the usual cumulant expansion 

(e» bAt ^ = l + bAtE 1 + ^E 2 (bAt) 2 + ^E 3 (bAt) 3 + --- 

= exp^A^ + ^^(bAt) 2 + Es ~ 3E f 2 + Ef (bAtf + ■ ■ .) (177) 

Here, b stands for either Va or —aa. Thus, for instance the average over ji + assumes the 
form 

[j/expL iaii ya ai! Am = ex pfe^ AtE +1 Va aii + X - (E +2 - E 2 +1 ) (Vh^At) 2 J 

E+3 — 3E + iE +2 + E\ 



X 



ex p EE 



6 



' +1 (Va aii At) 3 + 



(178) 



In the continuum limit the function a a i is replaced by the field ip + (t, x), whereas in the limit 
V — > the expression a a ^/V gives rise to the field ip(t,x). The sum over a together with 
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At gives rise to the time integral and the sum over i together with the volume element leads 
to the spatial integral £V V ~ > I dx. In the first term of the exponential in (178) we thus 
obtain 

^2^2AtE +1 Va a>i ^E +1 f dt J dxi) + {t,x) . (179) 

a i 

The continuum limit for the cumulants of second and higher order is not so obvious. We 
assume the simplest nontrivial distribution for /j,±, in which only the variance term has a 
finite limit, when At — > and V — > 0, whereas the contributions of higher-order cumulants 
vanish, for instance 

(E +2 - E 2 +1 )VAt ^ a+, At^O, V^O, (180) 
(E +3 - 3E +1 E +2 + Eli) (VAtf ->■ , At ->• , V ->■ . (181) 

Therefore, the contribution of the average over /i + to the effective dynamic action assumes 
the form 

S + = jdt jdx ^E +1 ^ + (t,x) + [^ f (t,a;)] 2 | . (182) 
For the average over /i_ a similar argument yields 



S_ = Jdt Jdx S^-E_ 1 ip\t,x)tp(t,x) + ^a_ [^{t, x)ip(t, x)] ' 



(183) 



These contributions to the effective dynamic action may, of course, be generated by suitably 
chosen normal distributions of /i±. 

This way of introduction of random sources and sinks has the annoying feature that it does 
not conserve the number of particles in the system. For a comparison with the treatment of 
this problem in the Langevin approach the random sources and sinks should be introduced 
in such a way that the particle number is conserved. The simplest way how to deal with this 
problem is to add to the random source a term proportional to the particle number, i.e. to 
use the "reaction constant" fi + V + /ii + n instead of ji + V in the master equation. The source 
terms on the right-hand side of the master equation (169) in this case assume the form 

= n + V[P(t,n-l)-P(t,n)} 

+ n 1+ [{n - l)P(t, n - 1) - nP(t, n)] . . . (184) 
The new part of the master equation corresponds to a branching process [4]. 
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The added term gives rise to the following contribution to the Hamilton operator 

H g2 (a + ,a) = -fi 1+ (a + -I) a + a. (185) 

Performing the steps described above we arrive at the contribution to the dynamic action in 
the following form 

S 1+ = Jdtjdx !^E 1+1 ^ + 1) ip + \°i+^ 2 (^ + 1) V} • (186) 

Now it is easy to see, that if we exclude the plain source (i.e. letting E +1 — a + — 0) and 
choose E\ + \ = E-i, the empty state remains absorbing one and the "mass term" oc ip^ip 
disappears in the dynamic action and we arrive at the dynamic action of random sources 
and sinks 

S gc = Jdtjdx jiW t2 V> + \o- (^V) 2 + \<ri+^ 2 + l) V}, (187) 

which conserves the average number of particles. 

The effects of the high-order terms are drastically different in the two cases amenable for 
a scaling analysis with the aid of the renormalization group. The time derivative term in 
the dynamic action 

S = - Jdt Jdx i)\t, x)d t ip(t, x) + ... (188) 

must be dimensionless in order to have nontrivial dynamics. Therefore the total scaling 
dimension of the number-density operator if)^(t, x)ip(t, x) is equal to the dimension of space 
and thus is positive. 

First, if the scaling dimension of the field ^ is equal to zero, d^ = 0, then the dimension 
of the field ip is positive (more precisely d^ = d) and the operator monomials in the second 
and third terms in (187) have the same scaling dimension. Since they are carrying the factor 
■?/> 2 , their scaling dimension is larger than that of ip^ip. Therefore, the second and third 
terms in (187) are IR irrelevant and should be discarded in the asymptotic analysis. 
Second, if the scaling dimensions of both fields are positive, then in the operator monomials in 
the second and third terms in (187) there is at least one "excessive" field factor in comparison 
with the first term, which renders them irrelevant. Thus, in these cases the IR relevant 
dynamic action of random sources and sinks reduces to the single term 

S' 9 c = J dt JdxE 1+1 ^ 2 ij, d^=0 V d i>t >0, d i ,>0. (189) 
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Third, if the scaling dimension of the field ip is zero, the scaling dimension of the field ift is 
positive and terms with "excessive" powers of ^ are IR irrelevant. So the starting point for 
the subsequent RG analysis is the source and sink action in the form 

SI = Jdtjdx \e 1+1 ^ 2 ^ + l - (<j_ + (V^V) 2 } , <h = . (190) 

3.2. ANNIHILATION PROCESS WITH RANDOM SOURCES AND SINKS 

Let us analyze the dynamic action of the diffusion-limited annihilation reaction A+A — > 

S 1 = -Jdtjdx ^d t ij - D ^V 2 ij + A A) W + (V^) 2 ] ^ 2 }+™o J dx^ f (x, 0), 

(191) 

from the point of view of scaling behaviour sketched in section 3.3.1. 

In the first case with d^+ = the nonlinear terms in action (191) are of equal scaling 
dimension. However, the source-sink part (189) is linear in the field ip with positive scal- 
ing dimension in contrast to the quadratic ip terms of (191). Therefore, the IR relevant 
interaction above two dimensions is (189) and the corresponding dynamic action is 

This dynamic action does not bring about any graphs with closed loops of the density prop- 
agator, which implies suppresion of the fluctuation effects. However, the scaling dimension 
of the interaction term is negative and may compensate for the positive dimensions of the 
irrelevant interaction terms. Therefore, the rest of the interaction terms are in fact danger- 
ous irrelevant operators and in this case definitive conclusion about the IR relevant action 
cannot be reached on the basis of the analysis of the scaling dimensions. 

In the second case with d^t > and d^ > the fourth-order term in action (191) becomes 
irrelevant. Either of the remaining third-order terms alone does not generate loops, therefore 
density fluctuation effects are brought about only, when both fields have the same scaling 
dimension d^t — d^ — d/2. In this case the IR relevant dynamic action is 

S IR2 = -Jdtjdx S^P + d t ^ - D ^ + V 2 ^ + 2A A)^V - £i+i^ + V J+no J dx^ + (x, 0). 

(193) 

Here, the scaling dimension of both interaction terms is (d/2) — 2 and vanishes at the critical 
dimension d c = 4, at which the dimensions of all the other interaction terms are positive 
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and they are unambiguously irrelevant. Effective action (193) is the dynamic action of the 
Gribov process [40], also known as the Reggeon model. Effects of random drift in this case 
with the use of the Obukhov-Kraichnan compressible velocity field have been analyzed in 
[41]. 

In the third case with = the fourth-order term in action (191) becomes irrelevant as 
well due to the positive dimension of the field ip + . By the same token, however, both terms 
of the source-sink action (190) are also irrelevant and we arrive at the IR relevant dynamic 
action 



SlR3 — 



dt J rfx ^dtiJ - D ^V 2 ij + 2AoA)^ 2 }+™o J rfx^ + (x, 0) . (194) 



An argument similar to that used for (192) shows, that the scaling analysis with this choice 
of field dimensions does not allow to resolve relevance of interaction terms. It should be 
recalled that the scaling dimensions of auxiliary quantities and the asymptotic behaviour of 
individual graphs is actually independent of the choice of the values of the field dimensions. 
Therefore, the effective action (193) with unambiguous classification of relevant and irrele- 
vant interaction terms describes the critical scaling behaviour amenable to the RG analysis. 

In summary, if the sources and sinks are chosen such that they conserve the mean number 
of particles in the system, the anomalous scaling behaviour in the system is that of the Gribov 
process. 

A different situation arises, if the plain source term is included into the analysis. Then 
there is a possibility that the system does not tend to the absorbing empty state but to an 
active state with a finite concentration of particles. In this case the starting point is the 
dynamic action with all the terms quoted above, i.e. 

S = Jdtjdx l-^dtip + D ^V 2 ij - \ D [2^ + (^) 2 ] i? + E +1 ^ 

+ \°+ (^) 2 + ^i+i^ + 1) ^ + ^i+^ +2 1 " + I)' ^ - E.^H 

+ ^<7_(V>V) 2 j+n Jdx V> f (x,0). (195) 

The stationarity equation brought about by this dynamic action for the field ip is (the 
stationary value ^ — as usual) 

dti) - AjV 2 ^ = -2A ^ 2 + E +1 + E l+l ^ - E.^p . (196) 
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However, the action expanded around the stationary value is rather complicated. To keep 
expressions simple, continue to consider the case E 1+1 = E_ x . Then the re-expanded action 
is 

S = JdtJ d^-^d t i) + D ^V 2 i) - V8y/E +1 \ D ^4> 

( ~E + i E_i^jE + i X D E + i<y l+ E + ia_ a + \ ,2 
^ 2 V2X D 4A 0J D 4A A) 2 1 ' 

E +1 cr 1+ ip^ E + ia l+ ^ A V2 ^E +1 X D a 1+ t/j^i/j 

2A A) 4A A) A 0J D 

It? ft Ft? — V~FT > a/^+i A £>o o~i+ a/ £+i A £>o C- \ .+2 



V E +i Aq-Pq ^i+ ^ ± , n x n ^ (J i+^ (J -\ ,t 2 /2 

+ g 1+ + ^ r }+noJ dx^(x, 0) . (197) 

In the critical limit — >■ 0. Since it is the expectation value of a nonnegative random 
quantity the variance <7+ vanishes as well. In the vicinity of the critical point we keep 
only the leading E + i and a + putting them equal zero in terms, where they are subleading. 
This simplifies the action a little bit 



s = JdtJ d*[-^d t i, + A)^vV - Vsy/E +1 \ D ^ + (^/== + y) ^ 

E +1 a 1+ ip^ E +1 a 1+ ip^ 4 ,2, V2 a/ E +1 X D a 1+ i/j^i/j 
+ hA-i^ IP H - 

y/E +1 X D a 1+ ^ip / 1 / 2 1 ( x n , a i+ , °"- \ /t 2 /2 

+ a 1+ + ai+ ^ r }+n J dx^(x, 0) . (198) 

Dimensional analysis of the canonical dimensions then yields the following cases. In the 
nonlinear parts without the critical parameters E + \ and a + the previous arguments hold, 
but in terms having powers of these parameters as coefficients the positive scaling dimensions 
of them must be taken into account. The free-field part of the action (198) suggests that the 
canonical dimension of E +x is four. In fact, the canonical dimension of a + remains a free 
parameter. 
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Proceeding in the same manner as above, we arrive at the following effective actions for 
the IR scaling limit. In the first case with d^t = the third and fourth powers of ^ and 
independent of ip or first order in tp are irrelevant (due to the coefficients proportional to 
E + \ or its square root) compared with terms oc ip^ 2 in action (198). Nonlinear in ip terms 
are irrelevant against the linear terms due to the positive dimension of i/j. Therefore, the IR 
effective action in this case is 

S = J dtj dx{-^d t ip + D ^V 2 tp -2V2 y/E +1 AoA^V 

+ + y) ^ + E. 1 ^}+n jd^^0). (199) 

Again, the interaction term remaining after the formal dimensional analysis does not bring 
about loops, although here we have a nontrivial correlation function of the field vp. The 
scaling dimension of this term is negative, however, rendering the irrelevant terms dangerous 
and prohibiting any definitive conclusion about the relevance of individual interaction terms. 

In the second case with d^t > and d^ > higher powers than the leading corrections to 
the free-field action of both fields are irrelevant. This argument leaves us with the dynamic 
action 

S = J dtj dx + D ^V 2 i) -2V2 y/E +1 A D ^V 

+ ( Ez^0^ + ?±>l ^ 2 + E_ 1 ^ 2 ij-2\ D ^ij 2 }+n JdxtffrO). (200) 

Contrary to the case discussed above, here the interaction term —2 AoA) ip^ip 2 generates 
loops alone due to the presence of the correlation function of the field ip. Therefore, two 
effective actions with nontrivial fluctuation contributions are possible. 

a) d^i > d^p. To keep the correlation function of the field ip for the loops, the variance er + 



must have a dimension less than that of y/E +1 . This yields the effective action 

S = J dtj dx^-^dtip + D ^V 2 tp - 2y/2y/E +1 A A^V 

+ ^±^ 2 -2X D ^^+n y"dx^+(x,0) (201) 

with the critical dimension depending on the scaling dimension of a + in the spirit of the 
description of tricritical scaling behaviour [19]. The model is logarithmic at six dimensions, 
however, because apart from the coefficient of the oc ip^ 2 the action is that of critical dynamics 
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of the ip 3 model. The upper critical dimension is determined by the scaling behaviour of a + 
in the critical limit a + —¥ 0, E +1 — > 0. 

b) — d^ — d/2. Both third-order terms are relevant and the effective action is basically 



(200). In this case the dimension of a + is larger than that of y/E + i and for simplicity we 
omit er + . Thus, the effective dynamic action may be written as 

s = Jdtjdx{ ~i) ] d t i) + Ay0 f v V -2V2 ^E +1 \ D ^ V 

+ E 'l ^ 2 + g_! ^ - 2 AqDq ^ V 2 | + n /dx^(x,0). (202) 

V2A -L'o J </ 

Note that this is a dynamic action describing the Gribov process with a random source 
independent of the active agent density. That the rate of change of the density due to the 
random sink is proportional to a power of density is a natural assumption. The assumption 
that the rate of change of the density due to the random source is proportional to a power 
of density is not natural. Therefore, the dynamic action (202) possibly predicts a critical 
behaviour of the Gribov process different from that discussed in the literature. 

In the third case with d^\ > and d^ = we arrive at the effective action (201). 

The analysis of scaling dimensions shows that we may actually lift most of the restrictions 
on the probability distribution of the transition rates of the type (180) and (181). Indeed, 
even if the higher order cumulants are finite, the scaling dimensions of corresponding terms 
in the dynamic action grow with the order of the cumulant with the exception of the case, 
when the transition rate is independent of the agent density. 

The scaling-dimension analysis of relevant and irrelevant interaction terms presented 
above appears somewhat formal. In particular, the arbitrariness of the scaling dimensions of 
the fields is an irritating detail. Therefore, it is instructive to repeat the analysis with the use 
of standard power counting. Consider first the case of particle-number conserving sources 
and sinks. The dynamic action is then (195) at the critical point, i.e. with E + i = a + = 
and E 1+1 = E_i. 

S = JdtJdK^-^ + d t ip + D^ + V 2 ^-\ D 2^j + +(ij + ) 2 ij 2 

+ E 1+1 ^ + V + ^i+^ +2 (^ + + l) 2 ^ 2 + ^-(^) 2 }+^o y"rfx^ + (x,0). (203) 
The divergence index of a one-irreducible graph is 

5=(d + 2)L-2I, (204) 
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where d is the dimension of space, L the number of loops and I the number of internal lines. 
The usual conditions relate the number of loops, lines, vertices, external field arguments 
and as well as the number of vertices V^-, in which the first subscript denotes the 
number of fields ip and the second the number of fields ip + in the action (203): 

L = I + 1 - (V 12 + V 21 + V 22 + V 23 + V 24 ) , 

I = V 12 + 2V 21 + 2V 22 + 2V 23 + 2V 2A - E,p , (205) 
/ = 2V 12 + V 21 + 2V 22 + 3V 23 + 4F 24 - . 

>From these equations it follows, in particular, 

E^ - E^ = -V 12 + V 21 - V 23 - 2V 24 . (206) 

Eliminating the number of lines / and the number of vertices V\ 2 we obtain 

5 = d + 2+(d- 4)V 21 + (d- 2)V 22 + dV 23 + (d + 2)V 24 - (d - 2)E i , - 2E^+ . (207) 

Multiplying relation (206) by an arbitrary coefficient a and combining with expression (207) 
we arrive at the following representation of the divergence index 

5 = d + 2 - aV 12 + (d - 4 + a)V 21 + (d - 2)V 22 + (d - a)V 23 

+ {d + 2- 2a)V 24 -(d-2 + a)E^ - (2 - a)E^+ . (208) 

Here, it is immediately seen that the choice of the value of the parameter a is tantamount 
to choosing the values of the scaling dimensions of the fields. Since the divergence index 5 
is independent of a, we conclude that the choice of the scaling dimension of the fields has 
actually nothing to do with the UV or IR behaviour of a one-irreducible graph. There is no 
mass parameter in the model, therefore a divergence index 5 with at least one negative coef- 
ficient of a vertex number indicates potential IR divergence, while a positive 5 corresponds 
to the usual superficial UV divergence. 

For further discussion, let us denote the generic vertex as v nm if} n if} +m and refer to any 
particular interaction term by its coefficient function v nm . >From expression (207) it would 
appear that vertices v 23 and v 24 always make a positive contribution to 5 regardless of the 
space dimension, therefore they are definitely IR irrelevant. The vertex t> 2 2 is irrelevant above 
two dimensions, whereas the vertex v 24 gives rise to IR divergences below four dimensions. 
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Since there is no known regular way to cope with these directly, we have to rely on the usual 
connection between the IR and UV divergences in the logarithmic theory and extrapolate its 
results below the critical dimension with the aid of the e expansion. Thus, if the vertex v 2 \ 
is present, the best we can do is to carry out the RG analysis around the critical dimension 
d c = 4, at which the vertex v 22 is also irrelevant, and we are left with the Gribov process. 
>From the point of view of the previous scaling dimension analysis representation (207) 
corresponds to effective action (193) at the logarithmic dimension d = 4 with the field 
dimensions d^ — d — 2, d^,+ = 2 and dimensionless coupling constant of the vertex v i2 . 

On the other hand, with the use of relation (206) we may eliminate any one vertex number, 
except V22, which then leads to changes in the coefficients of the remaining vertex numbers 
(apart from V 22 ) and to a different classification of the relevance of a given vertex, although 
the index of any graph does not feel the change. For instance, choosing a = 4 — d in (208) 
to eliminate V 21 we arrive at 

5 = d + 2 + (d- 4) Via + (d- 2)V 22 + (2d - 4) V 23 + (3d - 6) V 24 - 2E^ - (d - 2)E^+ . (209) 

Here, v 22 , v 2 % and v 2 4 are all irrelevant at d > 2, whereas V\ 2 is relevant at d < 4. Thus, if 
V12 > 0, the only critical regime amenable to an RG analysis is that of the Gribov process. 
It is evident from the preceding discussion that eliminating a vertex number from the 
expression for the divergence index is tantamount to putting the scaling dimension of that 
vertex equal to zero and thus fixing the scaling dimensions of the fields correspondingly. The 
aim of the power-counting analysis is to discard all graphs of a perturbation expansion of 
a Green function with the scaling dimension (i.e. the divergence index) larger than that 
of the leading-order expression. The divergence index is independent of the choice of the 
scaling dimension of the fields. In fact, the only tunable parameter it depends on is the 
space dimension. Putting the scaling dimension of a vertex equal to zero means that we are 
looking for the critical behaviour of the model at a space dimension at which there are no 
restrictions on the number these vertices in the relevant and marginal graphs. 

Since we are analyzing the effect of random sources and sinks on the pair annihilation 
process, we should keep the v 2 \ vertex in the effective IR model at any rate and therefore, for 
the purposes of the physical model, put the scaling dimension of this vertex equal to zero. 
This means that d^ = 2 and d^,+ — d — 2 and the that the critical behaviour of the model is 
that of the Gribov process, which belongs to a universality class significantly different from 
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that of the pure annihilation reaction model. 

>From the point of view of classification of terms in power counting the two effective 
actions (192) and (194) of the previous scaling-dimensions analysis correspond to elimination 
of E^+ and E^, respectively, from the expression for the divergence index. This yields 

5 = d + 2 - 2V 12 + (d - 2)V 21 + (d- 2)V 22 + (d- 2)V 23 + (d- 2)V 2A - dE^ (210) 
= d + 2 + (d - 2)V 12 - 2V 21 + (d- 2)V 22 + (2d - 2)V 23 + (3d - 2)V 2A - dE^+ . (211) 

Relation (210) suggests that all other vertices than v\ 2 are irrelevant above two dimensions 
leading to effective action (192). However, all these are "dangerous" irrelevant operators 
in the sense that the negative dimension of the vertex V\ 2 may render the dimension of a 
graph with irrelevant operators negative anyway. Since all dangerous operators have equal 
dimensions, they cannot be classified by the degree of "danger ousness". Relation (211) leads 
to a similar phenomenon with respect vertices other than v 2 ±, which has negative dimension 
there. However, here irrelevant operators have different dimensions and this might serve as 
a basis for classification of some operators more dangerous than others. 

Let us recall that all expressions (207) - (211) for the divergence index are equal and 
reflect different choices of the ambiguous field dimensions. We see that in a multicoupling 
case it is rather difficult to arrive at the consistent conclusion about the relevance of different 
vertices solely on the basis of the analysis of scaling dimensions of fields and vertex operators. 
In particular, to put either of the field dimensions equal to zero at the outset appears to 
be quite misleading. Assuming both field dimensions nonvanishing, put them equal to each 
other to obtain 

5 = d+2+(l - 2^j V 12 +(^ - 2^j V 21 +(d-2)V 22 +(^ - 2^j V 23 +(2d-2)V 24 ~E^~E^ + . 

(212) 

Here, indeed, all coefficients of vertex numbers depend of the space dimension explicitly and 
immediately lead to the conclusion, that the critical scaling behaviour tractable within the 
RG is that of the Gribov process. 

Consider then the action (197) for the active state. The divergence index does not give 
the whole truth in this case due to the presence of the "temperature" parameter E +1 , but it 
reflects the contribution of the wave- vector and frequency integral anyway. Many new types 
of vertices appear, but the set of conditions imposed on their numbers remains the same. 
In this case it is convenient to analyze the IR and UV behaviour separately. The reason is 
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that the temperature parameter is not involved in the power counting of UV divergences in 
any other way that new vertices have appeared. Thus, the power counting goes in the same 
fashion as above. Denoting new vertices by tilde, we obtain 

L = / + 1 - (v 02 + V 03 + V 04 + V 12 + V13 + V u + V 2 i + V 22 + V 23 + V 24 ) , (213) 
I = V 12 + Via + V u + 2V 21 + 2V 22 + 2V 23 + 2V 24 - , (214) 
/ = 2V 02 + 3V 03 + 4V 04 + 2V 12 + 3V 13 + 4V 14 + V 21 + 2V 22 + 3V 23 + 4V 24 - E^+ . (215) 

Consequently 

Eip - E^+ = -2V 02 - 3V 03 - 4V 04 - Vi 2 - 2V 13 - 3V 14 + V 21 - V 23 - 2V 24 . (216) 
>From (213) and (214) it follows that 

6 = (d + 2)L - 21 = d + 2 - (d + 2) V 02 - (d + 2) V 03 - (d + 2) V 04 - 2V 12 - 2Vi 3 

- 2V 14 + (d - 2) V21 + (d-2) V 22 + (d-2) V 23 + (d - 2) V 24 -dEj,. (217) 

Adding the relation (216) multiplied by the coefficient a to expression (217) we arrive at the 
representation 

S = d + 2 + (2a - d - 2) V 02 + (3a - d - 2) V 03 + (4a - d - 2) V 04 

+ (a - 2) V12 + (2a - 2) V 13 + (3a -2)V u + (d-2- a) V 2l 
+ (d-2)V 22 + (d-2 + a) V 23 + (d - 2 + 2a) V 24 -a^-(d-a) E^+ . (218) 

Here, by the choice of the value of the parameter a we could try to find a representation 
convenient for the classification of the vertices. However, in this case the vertex with the 
negative dimension ?p + ip + , may be absorbed in the pair correlation function of the field ip 
included in the elements of the graphical representation. In calculation of the UV index 
the field dimensions are then fixed by the condition that the coupling constant of the term 
oc ip + ip + is dimensionless. There is no ambiguity in the expression for the index and denoting 
the number of correlation functions in a one-irreducible graph I we arrive at relations 

L = I + I + 1 - (v 03 + V 04 + V12 + Vi 3 + V u + V 2 i + V 22 + V 23 + V 24 ) , 
I + 21 = Vi 2 + Vi 3 + Vu + 2V 2i + 2V 22 + 2V 23 + 2V 24 - E^ , (219) 
/ = 3V 03 + 4V 04 + 2V 12 + 3V 13 + 4V 14 + V 21 + 2V 22 + 3V 23 + 4V 24 - E^+ . 
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>From here it follows that 
6 = (d + 2)L-2I-AI = d + 2+(^ + ljv m + (d + 2) V 04 

+ V ^ + dV ™ + ( y + x ) Vu + - 3 ) V * + (d- 2)V 22 

+ (y ~ l ) V23 + 2dV2A ~ (t ~ x ) ^ " (J + x ) • (220) 

The same relation follows from (218), when the parameter a is chosen such that the coefficient 
in front of V 02 vanishes. >From (220) we see that when the vertex v 2 i is marginal, all the 
rest are irrelevant and we arrive at the situation described by the action (201). 

This, however, is not the whole story, because in the case of IR behaviour we are interested 
in the limit of vanishing temperature parameter. Since coupling constants of new vertices 
are functions of E +1 , they affect the limit of small E +1 directly. The dependence of a graph 
on E + i in the critical limit (ut -> 0, k { — > 0, E +1 — >■ 0, U{ = O (y/E+i), kf = O (^E + i)) 
is readily estimated by a suitable scaling of variables in the loop integrals, if these integrals 
have negative dimensions, in which case the upper limit in integrals of a renormalized graph 
may be sent to infinity. The result of the scaling of integration variables is a power of the 
parameter E + \ multiplied by a function of reduced frequencies and wave- vectors 0Ui/y/E + i, 
kiE+l^ 4 . The wave-number integrals defining this function are UV and IR finite for all values 
of the reduced frequencies and wave- vectors and thus possess a finite limit, when the latter 
vanish. Therefore, any IR-singular behaviour is signalled by a negative overall power of the 
parameter E +1 , which takes into account both the scaling of integration variables (this gives 
the divergence index 5) and the additional powers of E +1 at the vertex factors of the graph. 
Thus, we might use the following "IR divergence index" brought about by the action (197) 
(here, a + is assumed to be subleading for simplicity) 

6m = 5 + 2V 02 + AVos + 4V 4 + 2V 13 + 2V 14 = d + 2 + ( % + 2 ) V 03 + {d + 2) V m 



- ( i - 2]v 12 + dV 13 + (^) V u +(%-2) V 21 + (d- 2)V 22 





(3d 


; + 




\J 




'3d 


+ 1 







+ (t " 2 ) V23 + {2d ~ 2)1/24 -\ E *-\ E ^- ( 221) 



to guide in the classification of the interaction terms. The same form may be obtained the 
generic expression (218) with the choice of a such that the coefficient in front of V 02 in S m 
vanishes. In order to keep the IR behaviour tractable, all coefficients of the vertex numbers 
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in (221) must be nonnegative. The lowest space dimension conforming to this requirement 
is the upper critical dimension of the model. We see from (221) that d c — 4 for model (197) 
and the effective action is indeed (202). 

CONCLUSIONS 

In conclusion, we have analyzed the effect of density and velocity fluctuations on the 
reaction kinetics of the single-species decay A + A — >■ universality class in the framework 
of field-theoretic renormalization group and calculated the scaling function and the decay 
exponent of the mean particle density for the four asymptotic patterns predicted by the RG. 

We have calculated the relevant renormalization constants at two-loop level and found 
the decay exponent of the mean particle density at this order of the e, A expansion for four 
IR stable fixed points of the RG, whose regions of stability cover the whole parametric space 
in the vicinity of the origin in the e, A plane. The decay exponent assumes the mean-field 
value in the basins of attraction of the trivial fixed point (125) and of the kinetic fixed point 
(133) with dominant fluctuations of the random force of the Navier-Stokes equation. At the 
kinetic fixed point with finite rate coefficient (130) the decay value of the decay exponent is 
determined exactly by the fixed-point equations. At the thermal (short-range) fixed point 
(127) the decay exponent possesses a non-trivial e, A expansion. We have calculated three 
first terms of this expansion. 

Using a variational approach, we have inferred a renormalized fluctuation-amended rate 
equation with the account of one-loop corrections. This non-linear integro-differential equa- 
tion has been solved iteratively in the framework of the e, A expansion and the scaling 
function for the mean particle density has been calculated for the four IR stable regimes. 
The scaling function assumes the mean-field form (exactly) in the basins of attraction of 
the trivial fixed point and the kinetic fixed point with dominant fluctuations of the random 
force. At the kinetic fixed point with finite rate coefficient and at the thermal fixed point 
the scaling function possesses a non-trivial e, A expansion, which we have calculated at the 
linear order. Fluctuations of the random advection field affect heavily the long-time asymp- 
totic behaviour of the system: the kinetic fixed points are brought about by the velocity 
fluctuations as well as the non-trivial series expansion of the decay exponent at the thermal 
fixed point (without velocity fluctuations, the decay exponent is fixed to the one-loop value, 
because there are no high-order corrections to the rate constant in this case). Predictions 
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of the renormalization-group analyses for the reaction A + A — > in quenched random 
fields have been corroborated by numerical simulations [10, 13]. In the case of dynamically 
generated random drift this is seems to be a much more demanding task, but would surely 
be highly desirable, since the experimental data for reaction processes is quite scarce. s 

We have investigated possible effects of random sources and sinks on the pair annihilation 
reaction A + A — > 0. Contrary to the frequently used approach, in which the sources and 
sinks are introduced into the Langevin equation, we have included them directly to the master 
equation, where their physical sense is clear. We have considered linear in particle number 
creation and annihilation reactions with random rate coefficients to model the sources and 
sinks. On the basis of the analysis of canonical scaling dimensions we have constructed 
effective actions, which are the starting point for an RG analysis of the critical behaviour 
of the systems under consideration. In all cases the effect of random sources and sinks to 
the large-scale, long-time behaviour of the Green functions is significant and changes the 
universality class of the model. Instead of the universality class of the pair annihilation 
reaction A + A — > , the asymptotic behaviour of the model with random sources and 
sinks belongs to the universality class of the Gribov process in the critical case and to a 
modified Gribov process in the critical limit of the noncritical model. In the former case it 
is demonstrated once again that the description of a stochastic process with the use of the 
Langevin equation is significantly different from the description in terms of a master equation. 
The random noise term in the Langevin equation corresponds rather to the account of effects 
of genuine random sources and sinks than to a description of the effect of microscopic degrees 
of freedom on the mesoscopic process. Here, the universality classes of the same reaction 
process are completely different in the case of the master equation without sources and sinks 
in comparison to the case of Langevin equation for the same process. 

In the noncritical case with a random source independent of the agent density the Gri- 
bov process is modified to account for effects in critical behaviour, when sources and sinks 
asymptotically vanish. The analysis of the dependence of scaling functions on the parame- 
ters of the probability distribution of sinks and sources in infrared limit is called for. This 
reminds the situation, which takes place in the theory of phase transitions, where statistical 
correlations of the order parameter depend on a "mass" (deviation of temperature from the 
critical value) and the dependence of the scaling functions on the '"mass" is investigated. 
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APPENDIX A: EXPLICIT FORM OF THE RENORMALIZATION CONSTANTS 



Lll 



A 



;i+0^ + (3e + 2) M + 6£ + l 
512m(1 + m) 3 £ "~ ; 
u + 5 



112 



(1 + Om 2 + (1 + 3Qm + 6£-4 
256m(1 + m) 3 (1 -0 



22 512m(1 + m) 3 ' 

B 11 = B 1 (u) + B 2 (u,Z), B 12 = A[B 1 (u) + B 3 (u,^], B 22 = -B^u) - B 2 (u, -1), 



where the functions B 1 , B 2 and B 3 are given as 

1 



B 2 (u,0 



-12m 3 (1 +u) In 



+ 32w 4 (l + w) 2 ln2 + 2w 3 x 



1024m 4 (1 + m) 3 (m- 1) |_ 1 + u 

(1 + m) 2 (m + 10) In - - 32m 3 (1 + w) 3 arctghi + 4(2m 6 + 6m 5 + 7m 4 + 10m 3 - 
3 2 



1 + 2m 



1 + 2u 



Au - 1) In .t ' "~ - 4m 3 (1 + m)[4m 2 + 20m + 9] In - — — + 16m 3 (1 + uf x 
,l + Mr 



2 + 2u 



u 



arctgh— — + 8m 3 (1 + uf(u - l)h + ^(3/2))m 2 (23m 4 + 38m 3 + 17m 2 + 

M + 1 



22m + 4) 



1 



128ttm(1 + m) 2 (m 
(2 + 4£)m 2 + (lOf + 8)m + 38 + 22£ 



dq J dzF(q,z,u) 



1024m(1 + m) 3 (2 + O 
(8f + 4)m 2 + (14f + 10)m + 22 - lOf 
1024m(1 + m) 3 
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with function F given by the expression 

F(X,Z,U) = (1 - Z ) 1 —. r- 

M(x,z,u) = (:t 6 + 1)[,2 3 (24m 3 + 24m 2 + 72m + 72)-,2(8m 3 + i2m 2 + 8m + 60)] + 

(x 5 + x)[-z 4 (40u 3 + 88m 2 + 120m + 264) + z 2 (-4u 3 + 16m 2 + 108m + 168) + 
4m 3 + 14m 2 + 28m + 18] + (x 4 + x 2 )[z 5 (16u 3 + 96m 2 + 48m + 288) + ^ 3 (12m 4 + 
64m 2 + 128m 2 + 96m + 180) - ^(4m 4 + 26m 3 + 92m 2 + 174m + 312)] + 
x 3 [-z 6 (32u 2 + 96) - ^ 4 (8m 4 + 64m 3 + 240m 2 + 144m + 600) + ^ 2 (-8m 4 + 
4m 3 + 84m 2 + 108m + 452) + 2m 4 + 6m 3 + 26m 2 + 58m + 36], 

N(x,z,u) = (1 + x 2 -2xz)(l+x 2 -xz)((l+u)x 2 + 2-2xz)(l + u + 2x 2 -2xz) 



C(m,0 = — C dz(l-z 2 )^G(z,u) + —^ (\ n ^^ + l + 

v ' w 8mtt J_ x v ; v ' ; 8m(1 + m) V 2m 



2 + m 



X 



u 



2m + 2 s 



(Al) 



where 



G(z,u) 



( 1 u)- ■ \uz' 2 | 2 
m(m + 3)z 

y/2u(l+u) ~ 



u — 1 2m 
In ■ 



2(1 + u)z 



2^2 



U"Z 



1 +M v/1^2 

ZU + M + 1 



7T 1 + Z 

arctan \ / 

2 V 1 - z 



+ 



n — arctan 



u) — u 



2^2 



arctan 



(2 + z)u 



^2u(l + u)-u 2 z 2 _ 

APPENDIX B: FIXED POINTS 



(A2) 



9*22(0 



8R 8192 



3VTf 3VT7 



BM) - 



432^17(1 + 2 (2 + 



(21384 - 648V17)r + 



(52512 - 2592V17)r + (22192 - 2736V17)f 2 + (72V17 - 29064)£ + 
720^ - 18768^ , 
64(i?(2 + 30-l) 



16(2 + 30 



27(1 + 243(1 + 4 (2 + 

64(1 + R) 16(2 + 30 
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Figure 1. The propagators of the model 




Figure 2. Interaction vertices describing velocity fluctuation and advection with corresponding 

vertex factors 
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Figure 3. Interaction vertices responsible for density fluctuations with corresponding vertex factor 




Figure 4. Two-loop graphs for the perturbation expansion of T^t^ 
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Figure 5. Two-loop graphs for the perturbation expansion of T^t^2 




Figure 6. Regions of stability 
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Table 1. Canonical dimensions for the parameters and the fields of the model 
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Table 2. Canonical dimensions for the (1PI) divergent Green functions of the model 
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